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ABSTRACT 


This  thesis  consists  of  an  interactive  program  that  enables  the  student  to  study  the 
orbital  motion  of  satellites  around  the  earth.  The  student  can  investigate  the  shape  of 
a  variety  of  orbits  by  varying  the  initial  position  and  velocity  of  the  satellite,  or  by  sup¬ 
plying  select  orbital  parameters  i.e.  initial  orbital  radius,  eccentricity,  and  inclination. 
Satellite  maneuvers  can  also  be  studied,  like  transfer  orbits  and  inclination  changes,  by 
command  velocity  changes  at  any  location  in  the  orbit.  Also  the  clTccts  of  the  perturb¬ 
ing  forces  due  to  the  oblateness  of  the  earth,  drag  for  low  earth  orbits,  and  gravitational 
attraction  from  the  sun  and  moon  can  be  investigated.  The  orbits  arc  displayed  in  either 
the  perifocal  coordinate  system  around  a  model  of  the  earth,  or  the  ground  track  can 
be  displayed  on  a  map  of  the  world.  Orbital  data  is  displayed  below  the  orbital  plot. 
The  display  is  enabled  by  the  use  of  display  integrated  software  system  and  plotting 
language  tDISSPLAl  subroutines. 


iii 


1  icoMsien  tn  ^  | 

ms  gram 

MIC  TAB 

Unatununcad 

Justlfioatl 

V ^ 

□ 

□ 

an _ .  -  „ 

By - 

Oistr 

Ibutloo/ 

|  Availability  Codas  { 

Dlat 

(vV 

Avail 

5p*« 

u 

0 

13 

TABLE  OF  CONTENTS 


I.  INTRODUCTION  . 1 

II.  PROGRAM  DESIGN  . 3 

III.  UNPERTURBED  ORBIT  . 7 

IV.  PERTURBED  ORBIT  . 10 

A.  ORBITAL  ELEMENTS  . 10 

B.  COMPUTE  PERTURBING  FORCES  . II 

1.  NON-SPHERICAL  EARTH  . II 

2.  ATMOSPHERIC  DRAG  . 13 

3.  PERTURBING  FORCE  DUE  TO  IIEAVENLV  BODY  . 14 

a.  SUN’S  POSIT. ON  . 16 

b.  MOON'S  POSITION  . lo 

C.  RATE-OF-CHANGi:  01*  ORBITAL  ELEMENTS  . 1? 

D.  NEW  ORBITAL  ELEMENTS  . 19 

V.  VELOCITY  CHANGES  . 20 


VI.  GRAPHICAL  PLOTS  . 22 

A.  PERIFOCAL  PLOT . 23 

B.  GROUND  TRACK  . 23 

C.  DATA  . 23 

VII.  CONCLUSIONS  AND  RECOMMENDATIONS  . 25 


APPENDIX  A.  ORBIT  PROGRAM 


26 


APPENDIX  B.  COORDINATE  SYSTEMS  . 77 

A.  'UK' :  GEOCENTRIC  -  EQUATORIAL  . 77 

B.  I'QW:  PERI  FOCAL  . 'S 


iv 


C.  'RSW :  ORBITAI . 79 

ID.  COORDINATE  TRANSI  ORMATIONS  . Sn 

APPENDIX  C.  ORBITAL  ELEMENTS  . S2 

APPENDIX  D.  SAMPLE  ORBITS  . S6 

LIST  OF  R' TERENCES  . 91 

INITIAL  DISTRIBUTION  LIST  . 92 


I.  INTRODUCTION 


A  visual  aid  for  students  new  to  orbital  mechanics  is  required  to  comprehend  fully 
the  dynamics  of  orbital  motion.  This  program  is  an  interactive  time  step  simulation 
program  that  calculates  and  plots  either  unperturbed  or  perturbed  elliptical  orbits.  The 
program  interacts  with  the  student  in  developing  the  initial  orbit.  Also  the  program 
enables  the  student  with  the  ability  to  change  the  velocity  of  the  satellite  at  a  specific 
location  in  the  orbit.  This  feature  will  permit  the  student  to  investigate  the  eifects  of 
commanded  velocity  changes  as  in  perigee  kicks,  apogee  kicks  and  inclination  chunges. 
The  user  can  also  modify  the  initial  position  and  velocity  of  the  satellite  at  the  com¬ 
pletion  of  any  orbit. 

The  student  is  given  an  opportunity  to  investigate  the  elTeas  of  perturbing  forces 
on  the  satellites  orbit  by  choosing  to  have  the  program  calculate  the  orbit  with  or  with¬ 
out  perturbing  forces.  The  variation  of  parameters  method,  as  seen  in  (Ref.  I:  pp. 

rj.  is  used  in  calculating  the  perturbing  orbit.  The  perturbing  forces  taken  into 
consideration  are  the  following: 

1 .  the  oWatcnos  of  the  earth 

2.  drag  for  low  earth  orbits 

3.  gravitational  force  of  the  moon 

A.  gravitational  force  of  the  sun 

In  order  to  review  fully  the  operation  of  the  program  (included  in  appendix  A)  and 
to  uncover  any  problems  or  limitations  that  plagued  the  programming,  the  program  has 
been  divided  up  as  follows: 

1.  program  design 

2.  unperturbed  orbit 

3.  perturbed  orbit 

4.  velocity  changes 

5.  graphical  plots 

The  programming  approach  and  equations  used  in  each  of  the  above  sections  will  be 
examined  in  there  respccti\e  chapters.  A  review  of  the  coordinate  s>  stems  used  and  their 
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transfcrmaiscr.s  hew  con  them  arc  included  in  appendix  B.  Since  all  the  equations  used 
m  the  calculation  oi  the  orbital  elements  are  from  reference  I.  they  will  not  be  reviewed 
in  each  chapter  but  will  be  included  in  appendix  C  for  a  quick  reference,  liquations  from 
other  sources  will  be  rcterenccJ  in  their  respective  chapters. 

Examples  of  perturbed  and  unperturbed  orbital  plots  for  a  variety  of  initial  orbital 
parameters  are  inc  v  J  in  appendix  D.  included  arc  plots  of  low  earth  orbits,  transfer 
orbits  and  geosynchronous  orbits. 


II.  PROGRAM  DESIGN 


In  designing  tins  program  an  attempt  was  made  10  make  it  noi  only  as  user  friendly 
as  possible,  bui  also  to  make  the  program  as  simple  as  possible  to  understand.  To 
achieve  these  coals,  the  program  would  have  to  be  written  in  a  logical  manner,  in  a 
computer  language  that  is  easy  to  follow,  the  program  would  have  to  run  on  terminals 
readily  available  to  students  (at  the  Naval  Postgraduate  Sehoo!  (NPSli,  and  the  program 
would  have  to  be  easily  used  by  students  with  a  minimum  amount  of  computer  or  orbital 
mechanics  knowledge. 

FORTRAN  was  chosen  as  the  programming  language  since  it  is  a  wildly  used  sci¬ 
entific  language  and  it  allows  for  very  structured  programming.  By  programming  in  a 
structured  -format,  the  program  can  be  expanded  in  the  future  with  a  minimum  amount 
of  time  required  to  understand  the  programming  code.  FORTRAN  also  allows  for 
double  precision  numbers  to  be  used  in  the  calculation  of  the  orbit.  This  is  critical  when 
round  olT  error  in  single  precision  could  be  greater  then  the  actual  change  that  one  is 
trying  to  model.  The  equations  in  the  descriptions  of  the  program  might  not  exactly 
match  tiie  equations  in  the  listings  because  of  special  programming  techniques  which 
must  be  included  in  most  computer  programs  to  handle  such  problems  as  "division  b\ 
/.ero". 

The  display  integrated  software  system  and  plotting  language  (DISSPLA )  package 
available  un  the  mainframe  computer  at  NPS  was  used  to  enable  a  variety  of  graphical 
displays  with  a  minimum  amount  of  programming.  DISSPLA  has  a  set  of  subroutines 
that  the  programmer  calls  to  display  data  contained  in  arrays.  This  requirement  forces 
the  program  to  load  arrays  with  the  satellites  position  in  order  for  it  to  be  plotted.  The 
TEC618  computer  terminal  and  associative  plotter  was  used  for  ease  of  gaining  hard 
copy  plots  of  the  orbits  and  the  diversity  of  locations  that,  arc  available  here  at  NPS. 
In  order  to  run  a  program  in  DISSPLA  the  user  must  first  define  storage  space  of  15«.«»k 
and  designate  temporary  disk  space,  and  then  call  DISSPLA  with  the  program  name. 
This  is  accomplished  with  the  following  commands: 

1.  DLITNC  STORAGE  15WK 
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i  cms 

?.  TDISK  4  DiS 
4,  DISSPLA  ORBIT 


To  nuke  the  program  user  friendly,  the  user  is  prompted  for  inputs  via  the  keyboard. 
The  entry  is  usually  a  number.  A  yes  or  no  response  can  be  entered  by  typing  "V"  or  a 
"N".  In  most  cases  the  program  does  a  check  to  sec  if  the  input  is  appropriate.  In  order 
to  make  u  as  easy  as  possible  for  the  student  to  get  the  desired  orbit  displayed,  the 
program  requires  only  the  initial  position  and  velocity  of  the  satellite.  The  initial  posi¬ 
tion  and  velocity  of  the  satellite  is  supplied  by  the  user  in  one  of  two  ways.  The  user  can 
input  the  position  and  velocity  of  the  satellite,  using  the  perifocal  coordinate  system 
(UK),  or  the  user  can  let  the  program  place  the  satellite  on  the  "1“  axis  of  the  UK  system 
at  the  radius  of  perigee  tRPj  distance  supplied  by  the  user.  This  latter  choiec  gives  the 
initial  location  of  the  satellite,  but  to  get  the  velocity  the  program  will  prompt  the  user 
for  one  of  the  following: 

1.  the  actual  velocity  in  the  IJK  system. 

2.  the  eccentricity  ic)  of  the  orbit.  In  which  case  the  velocity  is  calculated  from  the 
following  equations: 

/.*/> 

u  *  yirr  *  semi-major  axis 


r.XR  =  -  ~  =  enerev  mass 

2,; 


Where  n  *  MG 

M  *  mass  of  earth 
G  “  Universal  gravitational  constant 


v 


\ 


2i  EXR 


3.  the  radius  of  apogee  (RA)  The  velocity  is  calculated  by  first  calculating  the  ec¬ 
centricity  (e)  from' the  following: 

RA  -  RP 
C  “  RA  -r  RP 

With  the  eccentricity  the  same  equations  used  above  are  used  so  calculate  the  ve¬ 
locity. 


In  order  to  give  the  velocity  a  direction  the  inclination  (i)  of  tlv  orbit  is  required  from 
the  user.  The  following  equations  are  used  to  calculate  the  velocity  vector: 


V/ =  n.M 


» i  vi1*;;., 


V,  M  1  NJIUI* 

The  program  will  check  to  ensure  that  the  orbital  eccentricity  j$  less  than  1.0,  if  it  is  not 
then  the  program  will  reject  the  inputs.  After  the  initial  input  are  accepted,  the  program 
will  do  calculations  for  the  six  orbital  elements  required  to  describe  the  size,  shape  and 
orientation  of  the  orbit,  and  to  pinpoint  the  position  of  the  satellite  along  the  orbit  at  a 
particular  time.  This  classical  set  of  six  orbital  elements  are  as  follows: 

1-  a.  semi-major  axis. 

2.  c.  eccentricity, 
i.  inclination. 

A.  Cl.  longitude  of  the  ascending  node. 

5.  ci.  argument  of  perigee  passage. 

A.  T.ume  of  perigee  passage. 

The  program  actually  calculates  more  orbital  elements  than  the  six  classical  elements 
required  to  plot  the  orbit,  this  is  done  in  an  effort  to  make  the  program  as  robust  as 
possible.  This  will  add  in  the  ability  to  expand  the  program  in  the  future. 

If  the  satellite  is  not  Initially  at  the  perigee  point  then  the  satellite  r?  first  be 
steppcJ  around  to  the  perigee  point.  The  program  then  enters  a  loop  that  calculates  the 
orbit  from  the  perigee  point  through  one,  complete  orbit  around  the  earth  and  back  to 
the  perigee  point.  The  orbit  is  calculated  in  steps  of  2  times  pi  divided  by  an  integer,  i.c., 
2  times  pi  divided  by  oO.  This  step  size  was  used  to  ensure  a  smooth  orbit  for  display 
purposes  and  also  to  get  within  adequate  distance  to  the  perigee  point  or  other  location 
for  a  velocity  change.  After  the  loop  is  completed,  the  program  will  offer  the  user  a 
choice  of  the  following  plots  to  check  the  orbit: 

1.  perifocal 

2.  grounutrack 

The  program  then  goes  into  a  loop  offering  the  user  the  following  choices  until  the  user 
decides  to  end  the  program: 

I.  plot  another  view  of  the  same  orbit. 
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if  the  u<cr  wishes  to  plot  another  view  of  the  same  orbit  then  the  user  may  use 
this  choice  to  reenter  the  display  portion  of  the  program. 

2.  plot  the  next  orbit « perturbed  or  unperturbed!. 

I  o  plot  the  next  orbit  the  satellite  is  stepped  around  the  complete  orbit  either 
with  or  without  perturbing  forces  eilecting  the  satellite. 

3.  change  the  initial  conditions. 

The  program  goes  to  the  beginning  of  the  program  and  allows  the  user  to  change 
the  initial  position  and  velocity  of  the  satellite. 

A.  change  the  velocity  at  a  specific  location 

Step  the  satellite’around  to  a  specific  true  anomaly  and  make  a  velocity  change 
at  that  location. 

5.  clear  the  previous  orbits  from  the  plot. 

Clear  the  memory  of  all  the  previous  orbits  and  only  retain  the  current  location 
and  velocity  as  the  initial  position  and  velocity. 

Before  each  new  orbit,  the  orbital  elements  arc  recalculated. 


There  arc  several  common  assumptions  and  constants  used  throughout  the  program 
i.e.  all  bodies  are  considered  to  be  spherically  symmetric  (this  allows  these  bodies  to  be 
treated  as  though  their  masses  arc  concentrated  at  their  centers  (point  masses)),  other 
assumptions  will  be  covered  in  their  respective  chapters. 
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III.  l\NPERTl'RBED  ORBIT 


The  subroutines  that  calculate  the  unperturbed  orbit  arc  the  most  widely  used  sub¬ 
routines  in  the  entire  program.  These  subroutines  are  called  to  step  the  satellite  around 
to  the  perigee  point  from  the  user  supplied  initial  position  and  velocity,  to  calculate  the 
next  unperturbed  orbit,  and  for  any  \eloeity  change.  No  matter  which  of  these  sources 
supply  the  initial  position  and  velocity  the  program  calculates  the  unperturbed  orbit  in 
the  same  manner.  The  only  difference  is  where  in  the  orbit  the  satellite  is  initially  when 
these  subroutines  arc  called.  Before  the  unperturbed  subroutines  arc  called,  the  orbital 
elements  are  calculated. 

The  unperturbed  subroutines  arc  called  by  a  single  subroutine  "UNPRET  which  has 
the  following  basic  algorithm: 

1.  Increment  time  by  the  time  step  size  (DTI.  The  time  step  was  chosen  as  the  period 
divided  In  tiny  to  give  a  smoo'h  plot,  but  more  importantly  to  ensure  that  the 
satellite  is  within  an  acceptable  distance  from  a  specific  location  lor  a  velocity 
change.  The  angular  error  caused  by  the  step  size  can  be  as  much  as  PI  50  from 
the  desired  point  for  a  circular  orbit  and  will  increase  for  more  eccentric  orbits. 
Hus  error  becomes  a  factor  when  the  user  is  making  velocity  changes,  and  therefore 
it  will  be  covered  in  that  chapter  in  further  detail. 

2.  Calculate  the  new  elements.  The  calculation  of  the  new  elements  is  the  heart  of  this 
algorithm.  The  size,  shape  and  orientation  of  the  orbit  remains  unchanged.  What 
is  required  is  the  position  of  the  satellite  along  the  orbit  as  a  function  of  time.  The 
problem  becomes  a  matter  to  solve  'the  Kepler  problem -predicting  the  future  po¬ 
sition  and  velocity  of  an  orbiting  object  as  a  function  of  some  known  initial  posi¬ 
tion  and  velocity  and  the  time  of  (light  (Ref.  1:  p.  1SIJ.  An  algorithm  using  these 
principles  will  follow: 

a.  A  time  step  (DT)  is  added  to  the  time  of  fliaht(TF),  time  of  flight  is  the  elapsed 
time  since  the  satellite  passed  the  perigee  point. 

7T«  TV-  DT 

b.  The  new  mean  anomaly  (MA)  is  calculated  from  the  new  time  of  flight,  and  the 
mean  motion  (MM). 

MA  -  MM  x  TF 

c.  With  the  new  mean  anomaly  the  new  eccentric  anomaly  (EA)  is  calculated. 
Because  the  solution  to  the  Kepler  problem  (d/.l »  EA  -ex  sin( £.-!))  is 
transcendental,  an  iterative  solution  based  on  the  Newton  method  of  root  find¬ 
ing  is  used.  The  root  in  question  is  a  solution  to  the  equation 
(MA  -  FA  ti’x  sini/f.-l)  =  0) .  This  algorithm  takes  the  form  of(Ref.  1:  p.  222): 

I )  MA.  -  FA.  -ex  sinfZT.7.) 


/:.» -  i:.i,  - 


i\u  -  m  l;. 

( 1  —ex  com/*'.  {*11 


Where  ihi<  equation  is  applied  initially  to  E.-U  ■  MA  and  then  reapplied 
until  the  dilTerenee  between  MA  and  MA.  becomes  small  enough  to  be  ig¬ 
nored. 

d.  The  new  true  anomaly  (v,)  is  calculated  from: 

cos~,(e-coM£-n) 

V°  "  ccosf £.•!)*"  I 

Calculate  the  new  position  and  velocity.  The  position  and  velocity  arc  calculated 
in  the  perifocal  coordinate  system  (PQW).  The  PQW  system  uses  the  orbit  as  its 
fundamental  plane  and  therefore  requires  only  two  coordinate  to  specify  the  satel¬ 
lite  s  position  and  velocity.  The  r.  coordinate  is  by  definition  always  equal  to  zero. 
The  position  of  the  satellite  is  calculated  as: 

.y„,  *■  r  cos  v 
yA  ■  r  sin  v 


The  velocity  of  the  satellite  is  calculated  as: 

^  y !  -  sm  v.-l 

v  -  .  ■  tc  4-  cos  vf ) 

\  / 

V.  *»  u 

A.  Store  position  and  elements  in  arrays  for  plotting,  in  order  for  the  program  to  plot 
the  orbit  the  radius,  true  anomaly,  inclination,  and  argument  of  perigee  must  be 
stored  in  arrays.  The  use  of  these  arrays  to  plot  the  orbit  will  be  explained  in 
chapter  t*. 

5.  The  process  is  repeated  until  the  satellite  is  at  the  perigee  point  and  the  true 
anomaly  is  two  pi. 

The  procedure  used  to  calculate  the  unperturbed  orbit  leave  very  little  to  be  modified 
by  a  programmer.  The  only  choices  that  had  to  be  made  concerned  step  size,  how  to  tell 
the  UNFRET  subroutine  that  the  perigee  point  had  been  reached,  and  a  value  of  ac¬ 
ceptable  error  for  newtons  method.  For  the  unperturbed  orbit,  the  step  size  just  had  to 
be  small  enough  to  produce  a  smooth  plot  of  the  orbit.  Two  indicators  for  perigee  were 
used,  one  was  that  the  true  anomaly  was  greater  than  6.21  radians  (two  pi  equals  6.2S 
radians)  and  the  time  from  the  previous  perigee  point  will  be  greater  then  the  period. 
The  two  indicators  were  logically  'mid'  together  to  ensure  the  perigee  point  was  reached. 
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I  he  disparity  hctuccr.  two  pi  and  r*.2l  radians  is  due  i?  il.c  error  produced  in  ih*  >.ti* 
ell.te  net  beginning  the  *rhit  at  exactly  the  perigee  point  ,md  the  step  <i/e  used  go 
around  the  orbit.  Use  acceptable  m.*c  of  error  t  er  newtons  method  was  set  at  I 
because  lor  an  unperturbed  orbit  this  would  be  the  major  contributor  to  any  error  in  the 
orbit  and  the  magnitude  of  this  error  would  be  acceptable.  However:  in  a  perturbed 
orbit  there  are  other  factors  contributing  to  determining  the  acceptable  error,  and  these 


will  be  discussed  in  the  next  chapter. 
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IV.  PERTURBED  ORBIT 


The  perturbed  orbit  uses  the  same  basic  routines  as  the  unperturbed  orbit  in  step¬ 
ping  the  satellite  around  the  earth  with  one  major  difference,  the  perturbing  forces 
produce  a  time  rate  of  change  of  the  orbital  elements  that  must  be  applied  at  ca*h  t;  ,« 
stem  The  variation  of  parameters  method  is  used  to  determine  this  influence  of  the 
perturbing  forces  on  the  orbital  elements.  The  analysis  is  simplified  by  using  the  orbital 
coordinate  system  ’RSW',  as  explained  in  appendix  B.  The  basic  algorithm  is  as  follows 
(Ref.  I:  P.  407]: 

1.  At  /  *  u  calculate  six  orbital  elements. 

2.  Compute  the  perturbing  forces  and  transform  it  at  i  -  r,  to  the  'RSW*  SYSTEM. 

3.  Compute  the  time  rate-of-change  of  the  elements. 

4.  Calculate  the  change  of  elements  for  one  time  step,  and  add  tiie  changes  to  the  old 
values  at  each  step  to  get  the  new  elements. 

5.  From  the  new  values  of  the  orbital  elements,  calculate  a  position  and  velocity. 

0.  Co  to  the  step  2  and  repeat  until  'he  final  time  is  reached. 

The  stops  in  the  algorithm  will  be  explained  in  the  following  sections: 

A.  ORBITAL  ELEMENTS 

The  standard  orbital  elements  a.  c.  i.  ii.  u  and  T  (or  M)  will  be  used,  where 
a  ■  semi-major  axis 
c  *  eccentricity 
i  **  inclination 

Cl  *  longitude  of  ascending  node 
<a  *  argument  of  perigee 
T  ■  time  of  perigee  passage 

(M,  *  mean  anomaly  at  epoch  -  M  -  n(t  -  /,))•  The  elements  arc  calculated  only 
at  the  beginning  of  the  orbit  from  the  initial  position  and  velocity  vectors.  The  elements 
arc  then  changed  continuously  throughout  the  orbit  by  adding  the  changes  due  to  the 
perturbing  forces.  For  the  perturbed  orbit,  the  satellite  will  always  begin  at  the  perigee 
point.  This  is  done  so  one  complete  orbit  is  from  perigee  point  to  perigee  point. 
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B.  COMPUTE  PERTURBING  FORCES 

The  variation  os*  parameters  method  requires  that  the  perturbing  forces  be  calculated  at 
each  step  hi  the  orbit.  In  order  to  do  this  a  model  of  each  perturbing  force  must  be 
developed.  The  following  perturbing  forces  where  used  in  calculating  the  total  perturb¬ 
ing  force  effecting  the  satellite: 

1.  oblatcness  of  the  earth 

2.  atmospheric  drag 

3.  gravitational  attraction  of  the  sun 

4.  gravitational  attraction  of  the  moon 

The  magnitude*  of  these  forces  have  an  enormous  range  of  values  and  are  dependent 
on  the  distance  the  satellite  is  from  the  perturbing  body.  Figure  i  on  page  12  shows  a 
graphical  representation  of  the  magnitude  of  the  perturbing  forces  in  a  log-log  plot  of 
perturbing  forces  per  unit  mass  (Ref.  2:  p.  IV-6I).  The  model  of  each  of  these  forces 
follows: 

I.  NON-SPHERICAL  EARTH 

The  earth  is  not  perfectly  spherical,  but  bulges  around  the  equator.  The  polar 
and  equatorial  diameters  arc  12713.0  Km  and  12756.3  Km.  respectively,  the  oblatcness 
results  in  a  perturbing  force  per  unit  mass  with  these  components  in  the  RSW*  coordi¬ 
nate  system !  Ref.  3:  p.  M|: 


„  (  >  .2 
F,  - r —  ( 1  -  .?  sin'U)  sm  (m.,1) 


„  ( -3«/-r;)  .  ,  . 

F,  - t5 —  ( sm'(/)  sm(i/a)  cos (u,)) 

r 


( -Wri) 


( sin  (/)  cos(/)  sinCt/y)) 


The  variable  and  constants  of  these  equations  arc  defined  below: 

1.  Variables: 

a.  it,  -  the  argument  of  latitude  and  is  equal  to  the  true  anomaly  v,  plus  the  ar¬ 
gument  of  perigee  u. 

i',»  ~  -  w 

b.  r  =  the  radius  from  the  center  of  the  earth  to  the  satellite. 
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Figure  1.  Comparison  of  perturbation  magnitudes. 
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2,  Oiwaniv 

«s.  u  ■  t!:c  gravitational  parameter  of  the  earth. 

I  Aim1  I 

jt  -  34St*)1.2 — 


b.  J:  ■  ihc  second  harmonic  of  oblatcncss  coefficient.  determined  by  experimental 
observations, 

-  l.n$23£- 3 


c.  r,  «  the  mean  radius  of  the  earth, 
rt  *  6.37S2£3A’m 

2.  ATMOSPHERIC  DRAG 


The  formulation  of  atmospheric  drag  equations  are  plagued  with  uncertainties 
of  atmospheric  fluctuations,  frontal  areas  of  orbiting  object  (if  not  constant),  the  drag 
coefficient.  and  other  parameters.  A  fairly  simple  formulation  will  be  given  hero.  Drag, 
by  definition,  will  be  opposite  to  the  velocity  of  the  vehicle  relative  to  the  atmosphere. 
Thus,  the  perturbing  foree  is 


r 


2<n 


)  •  CD  ».IR»  DKS  >v*v 


The  velocity  vector  is  in  the  ‘1JK’  system  so  the  resulting  force  is  also  in  the  TJK*  sys« 
tern.  Therefore  a  transformation  to  the  'R$\V*  system  is  required. 

The  variables  and  constants  of  this  equation  arc  defined  below: 

1.  Variables: 

a.  v  *  speed  of  vehicle. 

b.  CD  ■  the  dimensionless  drag  coefficient.  The  drag  coefficient  CD  has  a  value 
between  l  and  2.  It  takes  a  value  near  1  when  the  mean  free  path  of  the  at¬ 
mospheric  molecules  is  small  compared  with  the  satellite  size,  and  takes  a  value 
close  to  2  when  the  mean  free  path  is  large  compared  with  the  size  of  the  satel¬ 
lite.  The  drag  coefficient  will  be  modeled  with  CD  -  2  when  the  satellites  alti¬ 
tude  is  greater  than  550km  and  equal  to  1  otherwise.  (Ref.  4:  p.  295) 

c.  DEX  -  atmospheric  density  at  the  vehicle’s  altitude.  The  density  is  spherically 
symmetric,  and  will  be  modeled  using  exponential  steps  using  the’  parameters  in 
'I  able  1  on  page  14  and  the  following  formula  [Ref.  1:  pp.  423-424): 

dir)  -  . 
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Table  I.  ATMOSPHERIC  PAR  \METERS  AND  V.\LITS 


km*  :  .*  I  !\  |  .’ 

.h.-l 

».if*»  }  i  ::<i%«o 

ini 

i.222<im«2 

i 

15" 

l.«'I.*"3 

1 5"o/< 

55o 

3.0E-8 

550 

* 

|.nl54$4C-»7 

2.2lb9SE-o7 

I5U0 

mm 

■  I 

4|iki 

i«»r-i2 

2.  Constants  set  to  ivpical  values: 

a.  m  *  mass  of  the  satellite,  set  equal  to  100kg. 

h.  AR  "  the  sectional  area  of  the  vehicle  perpendicular  to  the  direction  of 
motion.  ** 

3.  PERI  )V£  TO  HEAVENLY  BODY 

The  satellite  v  ...„c  perturbation  force*  due  to  the  gravitational  clfccts 

of  the  sun  and  the  moon,  i  he  perturbation  force  from  a  perturbing  body  is  the  differ¬ 
ence  between  the  gravitational  force  due  to  the  perturbing  body  at  the  satellite  and  the 
graMtation.il  force  the  satellite  would  experience  if  it  were  at  the  center  of  the  earth. 
I  rom  l  igura  3  on  page  15.  the  perturbing  force  per  unit  mass  of  the  satellite  is 


The  variable  and  constants  arc  defined  below: 

1.  Variables: 

a.  rf  -  distance  from  the  earth  center  for  the  perturbing  body 

b.  if  -  unit  vector  from  the  earth  to  the  perturbing  body 

c.  r  *  distance  from  earth  center  to  the  satellite 

d.  i,  “  unit  vector  from  the  earth  to  the  satellite 

2.  Constants: 

a.  *  gravitational  constant  of  the  perturbing  body  -  MrG 

The  subscript  p  is  to  be  replaced  by  s  if  the  perturbing  body  is  the  sun.  and  by  m  if  the 
perturbing  body  is  the  moon.  We  will  assume  that  r  <  <  r,  then  the  equation  above 
becomes 
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•a 

The  unit  vectors  i,  and  if  can  be  written  in  terms  of  the  'IJK*  system  as: 
i,  m  ( cos(Q)  cos(t^)  -  sin(Q)  cos(i}  sin(tfj))/  +  (  cos(ty)  sin(Q)  +  cos(w)  cos {/)  sin(;^))J  + 

( sin  (i)  sinfifj))A* 


V  -  ( co$lQr)  cos(tf}p)  -  sin(flf)  cos(r^)  sin(^))7  + 

( cos(^)  sin(Qp)  -r  cos(Wp)  cos(/;)  sin(^))7  +  ( sin(^)  si n(ity))A* 

where  Q,  i,  and  u,  are  the  orbital  elements  of  the  satellites  and  Q„  i,,  and  %  are  the 
orbital  elements  of  the  perturbing  body.  The  formulas  above  use  the  TJK'  system,  and 
as  such  the  resultant  forces  must  be  transformed  to  the  'RSW'  system.  Models  of  the 

M* 

sun  and  moon  orbits  are  required  to  calculate  r,  and  rj.  The  models  used  in  the  program 
for  the  sun  and  moon's  orbits  follows:  [Ref.  3:  pp.  73-74) 
a.  SUN’S  POSITION 

In  order  to  model  the  suns  orbit,  a  number  of  simplifications  had  to  be 
made  in  the  actual  parameters  of  the  suns  orbit.  First  the  sun  will  be  assumed  to  be  in 
a  circular  orbit.This  means  that  the  radius  (r)  to  the  sun  will  be  constant,  and  the  ec¬ 
centricity  (c)  will  equal  0.0  instead  of  its  true  value  of  0.01 7.  The  other  assumption  will 
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he  10  place  the  sun  on  me  T  axis  of  the  IJK  system  at  the  beginning  of  the  program 
and  have  tt  progress  through  its  orbit  as  the  program  runs.  These  changes  will  not  clTcci 
the  perturbing  force  m  »i.y  noticeable  magnitude. 

The  following  variables  and  constants  where  used  in  the  program  to  model 
the  suns  orbit  after  applying  the  simplifications:  (Ref.  3:  pp.  75*7$) 

I.  Constants: 

C 

a.  Gravitational  Constant:  G  **  6.6r E  -  1 1  -V V- 

*«• 

b.  Sun's  Mass:  w#»  1 .99D0A.gr 

c.  Sun's  Gravitational  parameter: 


/i.-  1.32733£20-~- 

Aj? 


d.  Sun's  eccentricity:  a,  -  0.0 

e.  Radius  of  orbit,  assume  sun  is  in  circular  orbit:  r,  -  l.49£I  lm 

f.  Sun's  inclination:  si  ■  23.45  deg.  ■  4.092“9709d-0l  radians 

g.  Longitude  of  ascending  node:  C2,  -  «.u 

h.  Argument  of  perigee:  w,  »«M.l 
2.  Variables: 

a.  The  true  anomaly  of  the  sun's  position  as  a  function  of  the  time  the  satellite  has 
been  in  orbit: 

'  u<  n ^  "  350  x  24  x  3600  77 

Where  TT  *  true  time,  the  time  the  satellite  has  been  in  orbit  (sec) 

a*  *• 

b.  Sun's  Position  vector:  r  ■  r  cos  vitP  +  r  sin  \\Q 

-  V 

c.  I  nit  vector  from  the  earth  to  the  sun:  i,  -  — ~ 


b.  MOOXS  POSITION 

In  modeling  the  orbit  of  the  moon,  similar  assumptions  where  used  as  with 
the  sun.  The  moons  orbit  will  be  assumed  to  be  circular,  actually  the  eccentricity  is 
equal  to  0.055.  By  placing  the  moon  initially  on  the  T  axis  of  the  'IJK'  system  along 
with  the  sun.  the  gravitational  forces  of  the  two  bodies  will  combine  to  a  maximum. 
However;  since  the  moons  orbital  period  is  only  27.3  days,  the  moon  will  not  stay  in  this 
alignment  and  the  magnitude  of  the  combined  forces  will  vary  with  time.  The  inclination 
of  the  moons  orbit  is  not  constant,  but  drifts  between  1S.3  and  2S.6  degrees  in  ten  years. 
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Also  the  longitude  of  the  ascending  node  »£2i  oscillates  between  13  ar.J  *13  degrees.  To 
simplify  this  the  inclination  will  be  chosen  as  a  constant  23.5  degrees  and  the  longitude 
of  the  ascending  node  as  <ui  degrees.  For  the  time  period  involved  in  calculating  the 
perturbed  orbit,  these  assumptions  will  not  make  any  significant  difference. 

The  following  variables  and  constants  were  used  in  the  program  to  model 
the  moons  orbit,  after  applying  the  simplifications: 


I.  Constants: 


a.  Gravitational  Constant:  G  -  G.C»?£—  1 1 

b.  Moon's  Mass:  m,  -  "  35£22Ag 


(.Wi 


c.  Moon  s  Gravitational  Parameter:  ■  G.\L  ■  4.90C12 


d.  Moon's  eccentricity:  e„  -  0.0 

e.  Radius  of  orbit,  assume  moon  is  in  circular  orbit:  r„  -  3.$44£$Am 

f.  Moon's  inclination:  i  *  23.5deg.  «  4.ioi523"4fM  radians 

g.  Moon's  longitude  of  ascending  node:  £*  «  0.0 

h.  Moon's  argument  of  perigee:  u*  «  o.u 

i.  Moon  s  period:  T •»  2  *.3  days  (period) 

2.  Variables: 

a.  The  true  anomaly  of  the  moon's  position  as  a  function  of  the  time  the  satellite 

has  been  in  orbit:  v  JTH  -  7T 

2  >;  24  x  .’(hhj 

b.  Moon’s  position  Vector:  r  ■  r  cos  \y +  r  sin  ve,@ 


c.  Unit  vector  from  earth  to  moon:  L  -  — — - 

IM 

The  models  of  the  sun  and  moons  orbit  calculates  the  position  vector  in  the  TQW* 
system  and  therefore  the  position  vector  must  be  transformed  to  the  'UK'  system. 


C.  RATE-OFCHANGE  OF  ORBITAL  ELEMENTS 

The  derivations  and  equations  of  the  ratcs-of-changc  of  the  orbital  elements  arc 
contained  in  reference  1  pages  39S  to  40(5.  Therefore;  only  a  summary  of  the  actual  an¬ 
alytic  expressions  for  the  ratc-of-changc  of  the  parameters  in  terms  of  the  perturbations 
will  follow: 

1.  Ratc-ol -change  of  the  semi-major  axis: 


ri 


r.  <m  i  «>\  l  ~  t** 

— - 1  — ==  v.  -  c  — S; —  in 
K  \  1  “  ‘ 

Where  ’.•*  s%  the  mean  motton  of  the  satellites  erbit. 


\  a 

.  Raic-ef-ehangc  of  the  ceeemricity: 

*  .  \  l-*!  sm>„ _ v  ..  „ 

-r-t - n’l - K-C — — — 3C - ; - 'V, 

3.  Rate-of-ehange  of  the  inclination: 

A.  r  >c0<,/  y- 

,w  M\t\  I  -  e* 

4.  Rate-of-ehange  of  the  longitude  of  the  ascending  node: 


4uc. 

•//  L 


r  nm  Ufi 


n'ti\  1  -  e"  sin  i 

5.  Rate-of-ehange  of  the  argument  of  perigee: 
-41  -  ft/-—- 1,  h/-—  lt  *  it/- 7- 
Where. 


i  *  I  r 


-  v  1  -  «'*  COS  V 


(^.,-C^r]t.inv0(.+7-7!—  I]F, 


.  ,ht  .  r  -'•cousin II,  „„ 

('5r)*“c 

1  -  e* 

6.  Rate-of-ehange  of  the  eccentric  anomaly: 
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7.  Rate-of-chanyc  of  the  mean  anomaly: 
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reduces  10  the  following  for  vtrculur  and  ecliptic  orbits 
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Where  the  Rate-ci-ehange  of  the  mean  motion: 
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D.  NEW  ORBITAL  ELEMENTS 

The  change  of  each  element  is  calculated  by  multiplying  the  ratc-of-change  of  the 
element  by  the  time  step  (DTl.  The  change  in  the  orbital  elements  arc  then  added  to  the 
current  values  of  the  elements  to  give  the  new  orbital  elements.  With  the  new  elements 
calculated,  the  satellite  is  stepped  forward  and  the  new  position  and  velocity  are  calcu¬ 
lated  in  the  same  manner  as  the  unperturbed  orbit  (chapter  3).  Also  as  with  the  unper¬ 
turbed  orbit,  the  process  is  repeated  until  the  satellite  is  «t  the  perigee  point,  indicated 
by  the  time  of  flight  (TF)  equal  to  the  period  of  the  perturbed  orbit. 


V.  VELOCITY  CHANGES 


The  ability  of  the  student  to  change  the  velocity  of  the  satellite  at  any  position  in 
the  orbit  is  a  vital  element  in  this  program.  With  velocity  changes  the  student  can  in* 
vestigate  the  effects  of  varying  the  satellites  velocity  as  in  transfer  orbits  and  inclination 
changes.  In  order  to  simplify  the  program  the  unperturbed  orbit  is  used  throughout  this 
routine.  The  velocity  change  algorithm  used  in  the  program  follows: 

1.  Rotate  to  velocity  change  location. 

The  user  is  given  the  choice  of  changing  the  velocity  of  the  satellite  at  the 
perigee,  apogee  or  at  any  true  anomaly.  If  the  user  chooses  perigee  or  apogee  as 
the  change  locations,  the  true  anomaly  is  set  equal  to  zero  or  pi  radians  respect- 
full;. .  With  the  location  of  the  velocity  change,  the  satellite  is  first  stepped  around 
to  the  desired  true  anomaly.  The  stepping  is  identical  with  the  unperturbed  orbit 
with  the  exception  that  the  stepping  terminates  when  the  true  anomaly  is  greater 
or  equal  to  the  desired  true  anomaly.  With  a  step  size  of  one  fiftieth  of  the  period, 
the  satellite  is  actually  stepped  around  to  a  location  near  the  desired  location.  This 
variance  can  be  reduced  by  decreasing  the  step  size  but  this  would  increase  the 
computation  time.  This  error  will  be  a  major  factor  in  precise  calculations  of 
transfer  orbits,  or  any  other  orbital  maneuver  where  precise  velocity  changes  arc 
requireJ.  However:  this  program  is  not  a  tool  to  calculate  precise  orbital  maneu¬ 
vers  but  rather  a  learning  tool  for  the  student  to  get  a  >’’el  for  the  results  of  velocity 
changes  in  a  satellite’s  orbit. 

2.  Change  the  velocity. 

W*  ith  the  satellite  at  the  desired  location,  the  program  calculates  and  displays 
for  the  user  the  satellite’s  current  velocity,  escape  velocity  and  circular  velocity  (the 
velocity  required  to  circularize  the  orbit).  The  program  will  not  allow  velocities 
greater  than  or  equal  to  the  escape  velocity.  The  user  is  given  the  option  to  enter 
a  new  velocity  in  the  ‘UK’  system  or  to  change  the  magnitude  of  the  velocity  in  the 
orbital  plane.  If  the  user  chooses  to  change  the  velocity  in  the  orbital  plane,  the 
program  will  prompt  the  user  for  the  magnitude  of  the  velocity  change,  and  multi¬ 
ply  this  change  by  a  unit  vector  in  the  direction  of  the  satellites' velocity.  This  ve¬ 
locity  change  vector  is  then  added  to  the  satellites  velocity  vector,  to  calculate  the 
new  velocity  vector. 

3.  Calculate  new  elements. 

The  orbital  elements  are  calculated  with  the  new  velocity  vector  and  the  satel¬ 
lite’s  position  vector. 

4.  Complete  the  orbit. 

The  program  will  complete  the  orbit  to  the  new  perigee  point  using  me  satel¬ 
lite's  position,  new  velocity  and  new  elements.  There  arc  a  number  of  problems 
that  arise  if  the  satellite  is  just  stepped  around  to  the  perigee  point.  For  example, 
with  velocity  changes  in  the  orbital  plane  the  apogee  and  perigee  directions  can 
physically  swap.  This  is  a  problem  when  plotting  with  the  perifocal  coordinate 
system  because  the  .V,  axis  points  toward  perigee.  To  avoid  problems  like  this  the 
arrays  used  in  plotting  the  orbit  must  be  cleared  and  the  satellite's  current  position 
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ar.J  velocity  be  treated  a»  initial  conditions.  However,  to  compare  the  olJ  and  new 
orbits  there  i*  a  desire  to  retain  as  much  of  the  previous  orbit  as  possible.  The 
\ clouts  changes  where  divided  into  the  following  four  cases  to  handle  tiicse  prob* 
!c»k 

a.  C  hange  velocity  in  the  orbital  plane  at  the  perigee  point  with  the  new  velocity 
greater  than  the  circular  velocity.  The  perigee  point  will  remain  the  same  so  the 
satellite  is  stepped  around  using  the  unperturbed  subroutines. 

b.  Change  velocity  in  the  orbital  plane  at  the  perigee  point  with  the  new  velocity 
less  than  or  equal  to  the  circular  velocity.  The  perigee  and  apogee  directions 
will  switch  so  the  plotting  arrays  are  first  cleared  and  stored  with  the  current 
location  data.  Because  the  satellite  is  now  at  the  apogee  point  the  satellite  is 
stepped  around  to  the  perigee  point  storing  the  second  half  of  the  orbit.  The 
entire  next  orbit  is  calculated  and  stored  to  get  a  complete  orbit. 

c.  Change  velocity  in  the  orbital  plane  at  the  apogee  point  with  the  new  velocity 
less  than  the  circular  velocity.  The  perigee  and  apogee  directions  will  remain  the 
same,  so  the  satellite  is  stopped  around  to  the  perigee  point  completing  the  orbit. 

d.  This  last  case  catches  all  the  following  velocity  changes;  velocity  change  in  the 
orbital  plane  at  the  apogee  point  with  the  new  velocity  greater  than  or  cqual  to 
the  circular  velocity,  velocity  changes  at  any  other  true  anomaly  in  the  orbital 
plane,  and  any  velocity  change  out  of  the  orbital  plane.  The  plotting  arrays  arc 
cleared  and  stored  with  the  current  location  data.  No  matter  where  in  the  orbit 
the  satellite  is.  the  satellite  is  first  stepped  around  to  the  perigee  point,  and  to 
ensure  a  complete  orbit  is  plotted  the  entire  next  orbit  is  also  calculated  and 
stored. 
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VI.  GRAPHICAL  PLOTS 


The  program  provides  two  types  of  graphical  displays  of  the  orbit,  a  display  in  the 
perifocal  coordinate  system  and  a  display  of  the  satellite's  ground  track.  Each  display 
type  is  useful  in  observing  different  aspects  of  the  orbit.  The  perifocal  display  will  allow 
the  user  to  see  how  certain  orbital  parameters  change  with  different  initial  positions  and 
velocities,  and  also  how  the  parameters  change  with  velocity  changes  at  varying  posi¬ 
tions  in  the  orbit.  The  ground  track  will  enable  the  user  to  gain  an  appreciation  for  the 
physical  location  of  the  satellite  above  the  earth,  and  sec  how  the  orbital  parameter  af¬ 
fects  the  path  of  the  satellite.  The  ground  track  will  also  display  the  precession  of  a  se¬ 
quence  of  orbits.  Both  displays  plot  the  position  steps  to  give  the  user  an  understanding 
of  how  the  satellite  speeds  up  at  perigee  and  slows  down  around  apogee. 

The  DISSPLA  package  on  the  mainframe  computer  was  used  to  enable  the  plotting 
of  the  orbits.  The  versatility  of  plotting  subroutines  of  DISSPLA  makes  the  actual 
programming  of  the  orbit  a  simple  matter  of  initialling  DISSPLA  for  the  type  of  mon¬ 
itor  being  used,  setting  up  the  plotting  area,  initializing  the  axis  and  axis  scale,  and  then 
plotting  the  desired  curve  from  points  contained  in  arrays.  This  is  a  simplified  explana¬ 
tion  of  DISSPLA.  but  for  further  details  on  DISSPLA  programming  refer  to  the 
DISSPLA  user’s  manual  (Rel.  5).  DISSPLA  also  supplies  subroutines  to  draw  a  variety 
of  projections  of  the  world  and  (ill  the  projections  with  coast  lines,  latitude  lines  and 
longitude  lines.  Tnerc  are  a  couple  of  DISSPLA  requirements  that  did  require  special 
handling  in  the  program.  The  requirement  that  the  data  be  supplied  in  arrays  forced  the 
program  to  load  arrays  with  the  required  position  and  parameters  and  to  keep  a  counter 
for  the  number  in  the  arrays.  The  array  format  requires  the  size  of  the  array  be  specified 
in  the  beginning  of  the  program.  The  array  size  needs  to  be  large  enough  to  hold  a 
number  of  orbits,  but  not  so  large  as  to  waist  storage  space.  The  program  will  continue 
to  add  orbital  data  to  the  arrays  until  the  user  chooses  to  delete  the  previous  orbits.  If 
a  new  initial  position  and  velocity  is  entered  or  if  the  arrays  will  overflow  with  the  next 
orbit  the  arrays  will  automatically  delete  all  previous  orbits.  DISSPLA  also  requires  that 
all  data  be  in  single  precision  format.  The  program  calculates  all  orbits  in  double  pre¬ 
cision  in  order  to  limit  the  effect  of  round-off  error,  but  by  using  the  single  precision  data 
for  plotting  will  not  affect  the  accuracy  of  the  plot  in  any  way. 


Thc  subroutines  ',scd  to  display  the  orbits  will  be  eo\ered  in  the  following  three 
sections: 


A.  PERIFOCAL  PLOT 

Thc  plotting  of  thc  orbit  in  thc  perifocal  coordinate  system  is  thc  easier  of  the  two 
types  of  plots.  Since  thc  perifocal  coordinate  system  has  thc  orbital  plane  as  thc  fun¬ 
damental  plane,  the  only  requirements  to  describe  thc  orbit  in  thc  perifocal  coordinate 
system  arc  arrays  with  the  true  anomaly  and  the  radius  to  the  satellite.  To  give  thc  user 
a  sense  of  the  size  of  the  plot,  thc  axis  length  varies  with  the  eccentricity  and  semi-major 
axis  length.  Also  a  plot  of  the  earth  is  plotted  to  the  same  scale,  with  thc  pole  or  center 
of  thc  plot  on  thc  origin  of  the  axis.  Thc  latitude  of  thc  earth  at  the  center  of  the  plot 
will  vary  with  the  inclination  of  the  orbit.  This  plot  will  allow  thc  user  to  see  a  relative 


view  of  thc  satellite's  coverage  in  the  minus  'Z‘  axis  direction  of  thc  perifocal  coordinate 


system. 


B.  GROUND  TRACK 

Thc  ground  track  plot  is  a  very  complex  subroutine  compared  with  thc  perifocal 
plot.  Because  thc  ground  track  is  not  a  continuous  curve  a  procedure  to  handle  thc 
satellite  ending  at  one  end  of  thc  plot  and  wrapping  around  to  thc  other  end  was  devel¬ 
oped.  Thc  wrap  around  problem  is  avoided  in  most  orbits  by  plotting  the  orbit  in  seg¬ 
ments  with  the  following  two  rules.  Each  segment  begins  at  the  beginning  of  a  new  plot 
or  at  thc  edge  cf  the  plot  area,  and  ending  when  the  satellite  would  wrap  around  to  thc 
other  side  of  the  plot.  At  thc  beginning  of  a  segment  if  the  position  of  thc  satellite  is 
within  five  degrees  of  the  edge  of  thc  plot,  that  position  and  any  other  positions  within 
that  five  degree  boundary  will  not  be  plotted.  Thc  segment  will  end  when  the  satellite 
is  within  ten  degrees  of  thc  edge  of  the  plot.  Thc  above  restrictions  imposed  on  thc 
segments  of  thc  plot  will  not  substantially  affect  the  interpretation  or  usefulness  of  the 
plot.  The  ground  track  is  plotted  on  top  of  a  cylindrical  equidistant  projection  of  the 
world,  with  thc  world  coast  lines  and  a  longitude-latitude  grid  for  reference. 


C.  DATA 

Information  concerning  the  orbit  is  displayed  on  thc  lower  half  of  thc  plot.  Thc 
information  is  designed  to  supply  thc  user  with  enough  of  the  basic  orbital  elements  and 
other  parameters  affecting  the  orbit  to  be  able  to  evaluate  what  basic  type  of  orbit  thc 
satellite  is  in.  and  the  effects  of  velocity  changes  and  perturbing  forces  have  on  the  orbit. 
The  following  data  are  plotted:  inclination! i).  semi-major  axis  (a),  eccentricity  {eh  period 


23 


i  perl,  apogee  and  perigee  velocity  and  radius,  average  time  ratc-of*ehangc  of  orbital  cl* 
ement'.  and  the  average  magnitude  ol*  perturbing  forces  per  unit  moss. 


VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  program  supplies  the  student  with  an  interactive  tool  to  study  the  orbital  mo¬ 
tion  of  satellites  around  the  earth.  The  student  can  investigate  a  variety  of  orbits  by 
varying  the  orbital  parameters,  command  velocity  changes,  and  observe  the  effects  of 
perturbing  forces. 

The  student  is  provided  with  two  options  for  entering  the  initial  position  and  veloc¬ 
ity  of  the  satellite.  The  program  could  be  expanded  to  provide  the  student  with  the  ad¬ 
ditional  options  of  entering  cither  orbital  parameters  or  a  ground  observation  data  and 
have  the  program  calculate  the  initial  position  and  velocity  from  this  data.  Also  the 
student  is  limited  to  orbits  with  eccentricities  less  than  one  (elliptic  orbits).  The  program 
could  be  also  be  expanded  to  include  more  eccentric  orbit  for  Lunar,  interplanetary,  and 
missile  trajectories.  The  perturbing  orbit  is  calculated  for  orbits  around  the  earth  with 
relatively  small  perturbing  forces  in  relation  to  the  earths  gravitational  force.  This  fact 
will  cause  the  program  to  produce  false  results  if  the  student  tries  to  calculate  lunar 
trajectories.  Special  routines  would  have  to  be  employed  when  the  perturbing  force  (the 
moons  gravitational  attraction)  is  comparable  to  the  earths  gravitational  attraction. 
This  will  not  become  a  factor  for  studying  current  satellite  orbits  out  to  the 
geosynchronous  radius  of  4224 1.1  km. 

The  velocity  change  subroutines  move  the  satellite  to  a  location  close  to  the  desired 
location  before  a  velocity  change  is  imposed.  By  reducing  the  step  size  in  the  velocity 
change  subroutine,  this  error  could  be  reduced.  Precise  orbital  transfer  maneuvers  can 
be  modeled  by  reducing  this  error  caused  by  the  positioning  of  the  satellite  prior  to 
changing  the  velocity.  The  program  will  currently  provide  the  student  with  useful  plots 
for  gaining  experience  with  various  transfer  orbits  by  varying  the  magnitude  and  location 
of  the  velocity  changes. 

The  output  of  the  calculations  of  the  orbit  arc  arrays  loaded  with  the  satellite's  po¬ 
sition  and  select  orbital  parameters.  The  DISSPLA  subroutines  that  plot  the  points  arc 
not  unique.  The  program  would  become  portable  to  personal  computers  with  these 
graphics  subroutines  written  in  FORTRAN  and  included  in  the  program. 

A  final  recommendation  is  that  the  display  of  the  ground  track  could  be  modified 
to  show  ground  coverage,  number  of  satellites  in  a  constellation,  and  othor  elements 
necessary  for  planning  a  real-world  artificial  satellite  application. 


APPENDIX  A.  ORBIT  PROGRAM 

PROGRAM  CEBIT  0RB00010 

*  THIS  PROGRAM  IS  AN  INTERACTIVE  TIME  STEP  SIMULATION  OF  ORB00020 

*  SATELLITES  AROUND  THE  EARTH.  PERTURBED  AND  UNPERTURBED  ORBITS  ORB00030 

*  ARE  CALCULATED  AND  PLOTTED.  VELOCITY  CHANGES  ARE  ALSO  PERMITTED  0RB0Q040 


* 

AT  SPECIFIED  TRUE  ANOMALIES. 

ORB00050 

ORB00060 

it 

A  LIST  OF  VARIABLES  USED  BY  THE  MAIN  PROGRAM  FOLLOWS. 

ORB00070 

it 

A 

»  SEMI-MAJOR  AXIS 

ORBOOOAO 

it 

AL 

*  ARGUMENT  OF  LONGITUDE 

ORB00090 

it 

AP 

«  ARGUMENT  OF  PERIGEE 

ORBOOIOO 

* 

cm 

«  VELOCITY  CHANGE  LOCATION  TRUE  ANOMALY 

ORBOOIIO 

it 

DT 

»  TIME  STEP 

QRB00120 

it 

s 

•# 

*  ECCENTRICITY 

0RB00130 

* 

EA 

»  ECCENTRIC  ANOMALY 

0RB00140 

if 

El 

«  I  VECTOR  OF  ECCENTRICITY 

0UB00150 

it 

EJ 

«  J  VECTOR  OF  ECCENTRICITY 

0RB00160 

it 

EK 

■  K  VECTOR  OF  ECCENTRICITY 

ORB00170 

it 

FR 

«  R  VECTOR  OF  TOTAL  FORCE 

ORBOOIBO 

it 

FS 

*  S  VECTOR  OF  TOTAL  FORCE 

ORB00190 

it 

•M  » 

«  V  VECTOR  OF  TOTAL  FORCE 

0RB00200 

it 

H 

«  ANGULAR  MOMENTUM 

0RB00210 

it 

HI 

*  I  VECTOR  OF  ANGULAR  MOMENTUM 

0RB00220 

it 

KJ 

*  J  VECTOR  OF  ANGULAR  MOMENTUM 

ORB00230 

it 

HK 

*  K  VECTOR  OF  ANGULAR  MOMENTUM 

ORB00240 

it 

I 

*  INCLINATION 

ORB00250 
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IOPT1*  PERTURBED  OR  UNPERTURBED  OPTION 

ORB00260 

it 

iopt: 

OPTIONS:  PLOT  NEXT  ORBIT,  CHANGE  INITIAL  VALUES, 

0RB00270 

it 

CHANGE  VELOCITY,  PLOT  ANOTHER  VIEW  OF  ORBIT,  QUIT 

ORB002SO 

it 

LAN 

«  LONGITUDE  OF  ASCENDING  NODE 

0RB00290 

it 

LP 

=  LONGITUDE  OF  PERIGEE 

ORBO03OO 

it 

HA 

*  MEAN  ANOMALY 

0RB0031C 

it 

MM 

*  MEAN  MOTION 

0RB00320 

it 

MU 

«  GRAVITATIONAL  PARAMETER 

ORB00330 
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N 

*  ASCENDING  NODE 

0RB00340 

it 

NI 

*  I  VECTOR  OF  ASCENDING  NODE 

ORB00350 

it 

NJ 

=  J  VECTOR  OF  ASCENDING  NODE 

ORB00360 

it 

NK 

*  K  VECTOR  OF  ASCENDING  NODE 

0RB00370 

it 

NUM 

*  STEP  COUNTER 

0RB00380 

it 

P 

*  SEMI-LATUS  RECTUM 

0RB00390 

iV 

PER 

*  PERIOD  OF  ORBIT 

ORB00400 

it 

PI 

*  PI 

0RB00410 

it 

RA 

=  RADIUS  OF  APOGEE 

ORB00420 

i: 

RE 

■  RADIUS  OF  EARTH 

ORB00430 

it 

R 

=  OR3ITAL  RADIUS 

ORB00440 

it 

RI 

=  I  VECTOR  OF  ORBITAL  RADIUS 

ORB00450 

it 

RJ 

=  J  VECTOR  OF  ORBITAL  RADIUS 

ORB00460 

it 

RK 

=  K  VECTOR  OF  ORBITAL  RADIUS 

ORB00470 

it 

T 

=  TIME  COUNTER  IN  ORBIT 

ORB00480 

it 

TA 

=  TRUE  ANOMALY 

ORB00490 

it 

TDA 

=  TOTAL  CHANGE  IN  SEMI -MAJOR  AXIS 

ORB00500 

* 

TDAP 

=  TOTAL  CHANGE  IN  ARGUMENT  OF  PERIGEE 

ORB00510 
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TDE 

TDH 

TPI 

td:ia 

TDMM 

TDLAN 

T? 

TFDRA 

TFSA 

TFMO 

TFSU 

TL 


TOTAL  CHANGE  IN  ECCENTRICITY 

TOTAL  CHANGE  IN  ANGULAR  MOMENTUM 

TOTAL  CHANGE  IN  INCLINATION 

TOTAL  CHANGE  IN  MEAN  ANOMALY 

TOTAL  CHANGE  IN  MEAN  MOTION 

TOTAL  CHANGE  IN  LONGITUDE  OF  ASCENDING  NODS 

TIME  OF  FLIGHT 

TOTAL  FORCE  OF  DRAG 

TOTAL  FORCE  OF  EARTH'S  OBLATENESS 

TOTAL  FORCE  FROM  MOON 

TOTAL  FORCE  FROM  SUN 

TRUE  Long itudt  AT  EPOCH 

TRUE  TIME  SINCE  SATELLITE  HAS  BEEN  IN  ORBIT 

SATELLITE  VELOCITY 

I  VECTOR  OF  SATELLITE  VELOCITY 

J  VECTOR  OF  SATELLITE  VELOCITY 

K  VECTOR  OF  SATELLITE  VELOCITY 


A  LIST  OF  THE  ARRAYS  USED  FOLLOWS: 

AINRAY  »  INCLINATION 
APRAY  -  ARGUMENT  OF  PERIGEE 
RARAY  «  RADIUS 
RIP.AY  «  I  VECTOR  OF  RADIUS 
RJRAY  «  J  VECTOR  OF  RADIUS 
RKRAY  ■  K  VECTOR  OF  RADIUS 
TARAY  »  TRUE  ANOMALY 
TIMRAY  «  TIME 

A  LIST  OF  SUBROUTINES  CALLED  BY  THE  MAIN  PROGRAM  WILL  FOLLOW: 


CALCEL 

CHGVEL 

INPUTS 

INTSUM 

NEWELT 

KEVPOS 

NEWVEL 

OPTION 

PLOTS 

PRETUR 

STORE 

UNPRET 


CALCULATES  THE  ORBITAL  ELEMENTS 

ALLOW  THE  USER  TO  CHANGE  THE  VELOCITY  OF  THE  SATELLITE 

PROMPTS  USER  FOR  INITIAL  POSITION  AND  VELOCITY 

INITIALIZES  THE  SUMS  IN  THE  ARRAYS 

CALCULATE  NEW  ORBITAL  ELEMENTS  FROM  TIME  STEP 

CALCULATE  NEW  POSITION  VECTOR 

CALCULATE  NEW  VELOCITY  VECTOR 

GIVE  THE  USER  THE  OPTIONS  P«raitt«d  IN  THE  PROGRAM 

PLOTS  THE  ORBITS 

CALCULATES  THE  PERTURBED  ORBIT 

STORE  THE  POSITION  DATA  IN  ARRAYS 

CALCULATE  THE  UNPERTURBED  ORBIT 


BEGIN  MAIN  PROGRAM 

DOUBLE  PRECISION  PI.MU.RI.RJ.RK.R.VI.VJ-.VK.V.HI.HJ.KK^H, 

+  NI ,NJ,NK,N,P,EI ,EJ,EK,E,A, I ,LAN,AP>TA,AL,LP,TL,PER,EA, 

+  MM,MA,T,DT,TF,  FR,FS,FW,TT,CHTA,RA,VA,TEMPTA,RE 

DIMENSION  TARAY(SOO) ,RARAY(500) ,RIRAY(500) ,RJRAY(500) ,RKRAY(500) , 
+  AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* 1 , LOOP , YORN , ORLOOP 

PI  =  3. 1415926535S9794 


0RB00S20 

CRB00530 

CRBOOS&O 

0RB00550 

0RB00560 

0RB00570 

ORBOOSBO 

ORB00590 

0RBO0600 

ORB00610 

0RB00620 

0RB00630 

ORB00640 

0RB00650 

0RB00660 

0RB00670 

ORBOOSBO 

ORB00690 

0RB00700 

0RB00710 

0RB00720 

0RB00730 

0RB00740 

0RB00750 

0RB00760 

0RB00770 

QRB00780 

ORB00790 

ORBOOSOO 

0RB00810 

0RB00820 

ORB00830 

ORB00840 

ORB00850 

0RB00860 

ORB00870 

0RB00880 

0RB00890 

ORB00900 

0RB00910 

0RB00920 

OR900930 

0RB00940 

0HB00950 

ORBOO960 

0RB00970 

0RB00980 

0RB00990 

ORBOIOOO 

ORBQ1010 

ORB01020 

ORB01030 

ORB01040 

ORB01050 

0RB01060 

ORB01070 
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ML*  *  3.9S6012D+05 
RS  *  6. 3781^50+03 

USER  INTRO  TO  PROGRAM 
CALL  INTRO 


ENTERED  MAIN  PROGRAM  LOOP 
LOOP  ■  ' Y* 

IF  (LOOP  ,EQ.  '  Y* )  THEN 

INITIALIZE  STEP  COUNTER  AND  TRUE  TIME 
SUM  -  1 
TT  *  0. 0 

PROMPT  USES  FOR  INITIAL  POSITION  AND  VELOCITY 
CALL  INPUTS(RI,RJ,RK,R,VI,VJ,VK,V,MU,LOOP,PI) 

EXIT  PROGRAM 
IF  (LOOP  .EQ.  'N* )  THEN 
GOTO  10 
ENDIF 


CALCULATE  AND  STORE  ORBITAL  ELEMENTS 

CALL  CALCEL(RI,RJ,RK,R,VI,VJ,VX,V,EI,EJ,EK,E,A,I,LAS, 

+  LP,TA,P£R,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI ,HJ) 

CALL  STORE( RI , RJ , RK , R ,TA .RIRAY .RJRAY , RKRAY , RARAY ,TARAY , 

+  NUM,I,AP,AINRAY,APRAY,TT,TIMRAY) 


PRINT  DATE  FOR  USER  TO  REVIEW 


PRINT*,' VI 
PRINT* ,'VJ 
PRINT* ,'VK 
PRINT*,'  V 
PRINT*, 'RI 
PRINT*, 'RJ 
PRINT*, 'RK 
PRINT*,'  R 


VI, '  KM/S' 

VJ, '  KM/S' 

VK,  1  KM/S' 
V,‘  KM/S' 

RI, '  KM' 

RJ, '  KM' 

RK  '  KM' 
R,'  KM' 


PRINT*, 'ECCENTRICITY  «' ,E 
DEG I  *  SNGL(( 180.  O/PI )*I) 

PRINT* , ' INCLINATION  , LEG I , '  DEGREES ' 


PERHRS  «  SNGL( PER/3600.0) 

PRINT*,' PERIOD  *' ,FERHRS, '  HOURS' 

PRINT*, 'ARE  THESE  VALUES  CORRECT?  ENTER  "Y"  OR  "N"  : ' 
READ* , YORN 

CALL  EXCMS('CLRSCRN') 

IF  (.NOT.  YORN  .EQ.  lY')  THEN 
GOTO  20 


ENDIF 


CALCULATE  TIME  STEP  AND  SET  TIMER  TO  ONE  TIME  STEP 
DT  «  PER/50 
T  ■  DT 


STEP  SATELLITE  TO  PERIGEE  POINT  AND  RECORD 
IF  ((TA.  GT.  0.  063).  AND.  (TA.  LT.  6.  21))  THEN 
TT  =  TT  +  DT 


ORBOiOBO 
QRB01090 
0RB01100 
0RB01110 
0RB01120 
0RB01130 
0RB01140 
ORBOllSO 
0RB01160 
0RB01170 
ORBOllSO 
CRB01190 
0RB01200 
0RB01210 
0RB01220 
0RB01230 
0RB01240 
0RB01250 
0RB01260 
0RB01270 
ORB01280 
ORB01290 
0RB01300 
0RB01310 
0RB01320 
0RB01330 
0RB01340 
0RB01350 
ORB01360 
0RB01370 
0RB01380 
0RB01390 
ORB 014 00 
0RB014I0 
0RB01420 
0RB01430 
0RB01440 
ORB01450 
ORB01460 
ORB01470 
ORB01480 
ORB01490 
ORBOISOO 
0RB01510 
0RB01520 
0RB01530 
0RB01540 
0RB01550 
0RB01560 
0RB01570 
ORB01580 
0RB01590 
ORB01600 
0RB01610 
ORB01620 
0RB01630 


+  + 


CALL  NEVELT( MM , MA , E , EA , TA , TF , DT , ?I , PER ) 

CALL  N?GS(RI ,RJ , RK , R , LAN , A? , I ,  TA,A,E) 

CALL  NVEL(  E ,  ?  ,T.< ,  LAS , A? ,  I , VI , VJ , VK , V.MC) 

SUM  -  NCM  +  1 

CALL  ST08EC  81 , RJ , RK , 8 , TA , 8 1 RAY , RJRAY , RKRAY , RARAY ,7ARAY , 
NCM, I , AP , A INRAY , AP8AY ,  TT , TIMRAY) 

T  «  7  +  D7 
0070  50 

ENDIF 

CALCULATE  ELEMENTS  FROM  PERIGEE  POINT 

CALL  CALCEL(RI,RJ,RX,R,VI,VJ,VK,V,EI,EJ,EK,E,A,I,LAN, 

L?,TA,P£R,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

DT  ■  PER/50 
7  *  DT 


0RB0164O 

0RB01650 

0RB01660 

0RB01670 

0RB01680 

0RB01690 

CRB01700 

0RB01710 

ORBO1720 

0RB01730 

0RB01740 

ORB01750 

ORBOI760 

ORBOI770 

0RB01780 

0RB01790 

ORBOIBOO 

0RB01810 

ORB01820 


STORE  FIRST  Unparturbad  ORBIT  0RB01810 

CALL  UNPRET( DT , PER , AL , LAN , AP , I , R I , R J , RK , R , VI , VJ , VK , V ,  ORB01820 

MU,PI,H,A,E,N,TA,P,MM,MA,EA,TF,T,NCM,RIRAY,RJRAY,  0RB01830 
RKRAY,RARAY,TARAYfAINRAY,APRAY, TIMRAY, TT)  ORB01840 

0RB01850 

INITIALIZE  SUMS  FOR  FORCE  AND  ORBITAL  ELEMENT  CHANGES  TO  ZERO  ORB01860 
CALL  INTSUM(TFEA,TFSU,TFMO,TFDRA,TDI,TDA,TDE,TDMM,TDMA,TDLAN,  0RB01870 
TDH.TDAP)  0RB01880 

0RB01890 

PLOT  FIRST  UNPERTURBED  ORBIT  0RB01900 

CALL  PLOTS ( RIRAY , RJRAY , RKRAY , RAHAY ,TARAY , NUM , PI , I , LP , A , E ,TF ,  ORBO 1910 
AINRAY,APRAY, TIMRAY, TFEA,TFSU,TFMO,TFDRA, PER, TDI ,TDA,0RB01920 
TDE,TDMM,TDMA,TDLAN, TDH.TDAP, MM, MA, LAN, H,AP,R,V)  0RB01930 


BEGIN  NEW  ORBIT  OPTIONS 
IOrTl  »  1.  Unparturbtd  ORBIT 

■  2.  Ptrcurbad  ORBIT 
*  3.  QUIT 

I0PT2  »  1.  PLOT  NEXT  ORBIT 

■  2.  CHANGE  INITIAL  VALUES 

«  3.  CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  Anomaly 
»  4.  PLOT  ANOTHER  VIEW  OF  SAME  ORBIT 

ALSO  ASKED  IF  WANT  TO  CLEAR  ALL  PREVIOUS  ORBITS 

CALCULATE  ELEMENTS  AT  PERIGEE 

CALL  CALCEL(RI,RJ,RK,R,VI,VJ,VX,V,EI,EJ,EK,E,A,I,LAN, 

LP,TA,PER,EA,MA,AP,AL,TF,P,FI,MU,MM,N,H,HI ,HJ) 

CHECK  FOR  POSSIBLE  ARRAY  OVERFLOW 
IF  (NUM  .GT.  425)  THEN 

PRINT*, 'ARRAYS  ARE  FULL' 

PRINT*,’ PREVIOUS  ORBITS  WILL  BE  ERASED!’ 

NUM  =  1 

CALL  STORE(RI,RJ,RK,R,TA, RIRAY, RJRAY, RKRAY, RARAY, TARAY, 
NUM , I , AP , A INRAY , APRAY , TT , TIMRAY ) 

ENDIF 

PROMPT  USER  FOR  DESIRED  OPTION 


ORB01940 
ORB01950 
ORBO 1960 
0RB01970 
0RB01980 
ORB01990 
QRB02000 
ORB02010 
ORB02020 
ORB02030 
0RB02040 
0RB02050 
ORB02060 
ORB02070 
ORB02080 
ORB02090 
ORB02100 
ORB02110 
ORB02120 
ORB02130 
0RB02140 
ORB02150 
ORB02160 
ORB02170 
ORB02180 
OR502190 
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*■  +  +  +  + 


CALL  OPTIONC IOPTI, ICPT2 , SUM .RIRAY, RJRAY , RKRAY , RARAY , 
r  TAF.AY , AINRAY ,  APRAY .TIMRAY) 

Initial!.*:*  SIMS  FOR  FORCE  AND  ORBITAL  ELEMENT  CHANGES  TO  ZERO 
CALL  IN?SUM(TFEA,TFSU,TFHO,TFDRA,TOI , TDA , TOE , TDMM , TDMA , TDLAN , 

+  TOH.TDAP) 

SET  TIME  COUNTER  TO  ONE  TIME  STEP 
T  «  DT 

OPTION:  PLOT  THE  NEXT  ORBIT 
IF  (IOPT2  .EQ.  1)  THEN 

CALCULATE  AND  PLOT  UNPERTURBED  ORBIT 
IFCIOPT1  .EQ.  1)  THEN 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ, 

+  RK,R,VI ,VJ,VK,V,MU,PI ,H,A, 

+  E, N, TA, P, MM, MA,EA,TF,T, SUM, RIRAY, RJRAY, RKRAY, 

+  RARAY ,TARAY , AINRAY , APRAY .TIMRAY ,TT) 

CALL  ?LOTS( RIRAY , RJRAY , RKRAY .RARAY .TARAY , SUM , 

+  PI, I, LP, A, E, IF, AINRAY, APRAY, TIMRAY, 

+  TFEA,TFSU,TFMO,TFDRA,PER, 

+  TDI, TDA, TDE, TDMM, TDMA, TDLAN, TDK, TDAP, 

+  MM,MA,LAN,H,AP,R,V) 

CALCULATE  AND  PLOT  PERTURBED  ORBIT 
ELSEIFC IOPTI  .EQ.  2)  THEN 

CALL  PRETUR(DT,PER,AL,LAN,AP,I, 

RI,RJ,RK,R,VI,VJ,VK,V,FR,FS,FV. 

MU, PI ,H,A,E,N,TA,P,MM,MA,EA,TF,T,NUM, 

RIRAY, RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , 
TIMRAY, TT,TF£A,TF5U,TFM0,TFDRA, 

TDI , TDA , TDE , TDMM , TOMA , TDLAN ,TOH , TDAP) 

CALL  PLOTS( RIRAY , RJRAY , RKRAY , RARAY .TARAY, SUM , 

+  PI, I,LP, A, E,TF, AINRAY, APRAY, TIMRAY, 

+  TF£A,TFSU,TFMO,TFDRA,PFR, 

+  TDI, TDA, TOE, TDMM, TOMA, TDLAN, TOH, TDAP, 

+  MM,MA,LAN,H,AP,R,V) 

END  IF 

GOTO  THE  BEGINNING  OF  THE  PROGRAM  TO  CHANGE  THE  INITIAL  VALUES 
ELSEIF  ( IOPT2  .  EQ.  2)  THEN 
GOTO  20 

CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  ANOMALY  AND 
PLOT  THE  NEW  ORBIT 
ELSEIF  (I0PT2  .EQ.  3)  THEN 

CALL  CHGVEL(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R, 

+  VI,VJ,VK,V,MU,PI, 

+  H.A.E.N.TA.P.MM.MA.EA.TF.T.NUM.RIRAY, 

+  RJRAY, RKRAY, RARAY, TARAY, AINRAY, APRAY, 

+  TIMRAY , IT, El, EJ,EK,LP, HI, HJ, IOPTI, 

+  TFEA,TFSU,TFMO,TFDRA,  TDI,  TDA,  TDE,  TDMM, 

+  TDMA, TDLAN, TOH, TDAP) 

CALL  PLOTSC  RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 

+  PI.I.LP.A.E.TF, AINRAY , APRAY , TIMRAY , 


ORB02200 

0RB02210 

0RB02220 

ORB02230 

0RB02240 

ORB02250 

0X102260 

0X102270 

0X102260 

0X102290 

0X102300 

QRB02310 

0X102320 

0RB02330 

0X102340 

0X102350 

0X102360 

ORB02370 

ORB02380 

ORB02390 

0RB0240O 

0X102410 

0X102420 

0RB0243Q 

0RB02440 

0RB024S0 

0R192460 

0RB02470 

0RB02480 

0RB02490 

0RB02500 

ORB02510 

0RB02520 

0RB02530 

0RB02540 

0RB02550 

0RB02560 

ORB02570 

ORB02S80 

0RB02S90 

ORB02600 

ORB02610 

0RB02620 

ORB02630 

0RB02640 

ORB02650 

ORB02660 

ORB02670 

0RB02680 

0RB02690 

0RB02700 

0RB02710 

ORB02720 

0RB02730 

ORB02740 

ORB02750 
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+  TFEA ,7F5U ,TFMO,TFDRA  ,  PER , 

*  TDI.TDA.TDE.TDMM.TDMA.TDLAN.TDH.TDAP, 

*  MM ,  MA ,  LAN , K , AP , R , V) 

*  PLOT  ANOTHER  VIEW  CF  THE  SAME  CR3IT 
ELSEIF  (I0PT2  .EQ.  4)  THEN 

CALL  ?LOTS( RIRAY , RJRAY , RKRAY , RARAY ,TARAY .SUM , 

+  PI ,I,LP,A,E,TF,AINRAY,APRAY,TIMRAY, 

+  TFEA, TF5U,TFM0,TTDRA, PER, 

+  TDIJDA.TDE.TDMM.TDMA.TDLANjTDH.TBAP, 

+  MM,MA,LAN,H,AP) 

*  STOP  THE  PROGRAM 

ELSEIF  (IOPT2  .  EQ.  5)  THEN 
GOTO  90 

ELSE 

PRINT*, ' INVALID  ENTRY! ' 

GOTO  80 
ENDIF 

*  CHECK  IF  SATELLITE  Impacted  THE  EARTH  AND  GO  TO  THE  BEGINNING 
IF  (R  .  LE.  6450. 0)  THEN 

PRINT*, ’SATELLITE  WILL  IMPACT  THE  EARTH!!!  ' 

PRINT*, 'PROGRAM  WILL  RESET  TO  THE  BEGINNING' 

GOTO  20 

ENDIF 

*  GOTO  THE  TOP  OF  THE  OPTION  LOOP 
GOTO  80 

*  GIVE  THE  USER  A  CHANCE  TO  RECOVER  THE  PROGRAM 
90  PRINT*, 'THIS  IS  YCUR  LAST  CHANCE! ' 

PRINT'S' DO  YOU  WANT  TO  CONTINUE?' 

PRINT'S 'AND  GOTO  THE  Beiiinnin*  OF  THE  PROGRAM?' 

PRINT'S' ENTER  "Y"  OR  "N,f  :  ' 

READ*, LOOP 
PRINT’S  LOOP 
GOTO  10 
ENDIF 

*  DISSPLA  SUBROUTINE  TO  TELL  GRAPHICS  TERMINAL  PLOTTING 

*  SESSION  IS  DONE 
CALL  DONEPL 
STOP 

END 


SUBROUTINE  INTRO 

*  THIS  SUBROUTINE  WILL  GIVE  THE  USER  A  Brief  INTRODUCTION  OF  THE 

*  USES  OF  THE  PROGRAM 

PRINT'S  'THIS  PROGRAM  IS  A  GRAPHICS  DISPLAY  OF  Satellite  ORBITS.  ' 
PRINT*, 'YOU  WILL  BE  ASKED  TO  INPUT  THE  INITIAL  VELOCITY  AND' 
PRINT’S 'POSITION  VECTORS  OF  THE  Satellite.  THE  PROGRAM  WILL  ' 
PRINT’S  'THEN  CALCULATE  THE  ORBITAL  PARAMETERS  AND  THE  ' 


CRB02760 

0RB02770 

CRB027B0 

0RB02790 

GRB02800 

0RB02810 

0RB02820 

0X102630 

0XB02640 

0RB02850 

0RB02660 

ORB02870 

0X102860 

0X802890 

0X802900 

0X102910 

0X102920 

ORB02930 

0RB02940 

ORB02950 

0RB02960 

0RB02970 

0RB02980 

0RB02990 

ORB03000 

ORB03010 

0RB03020 

0RB03030 

0RB03040 

ORB03050 

0RB03060 

0RB03070 

0RB03080 

0RB03090 

ORB03100 

0RB03110 

ORB03120 

ORB03130 

0RB03140 

ORB03150 

0RB03160 

0RB33170 

0RB03180 

0RB03190 

ORB03200 

0RB03210 

0RB03220 

0RB03230 

OKB03240 

CRB03250 

0RB03260 

0RB03270 

ORB03280 

ORB03290 

ORB03300 

0RB03310 
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print*, 

v.iM'S 

PaINT*, 

PT-IST*, 
HINT'S 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT* , 
PAINT*, 
PRINT*, 
PRINT*, 
PRINT*, 
PRINT'S 
RETURN 
END 


’Unperturbed  CSBIT.  THE  USER  VILE  THEN  HAVE  THE* 

’CHOICE  OF  DISPLAYS: 1 

1  -FSRIFCCAL  ( SHOWS  RELATIVE  SIZE  OF  CR3IT)' 

’  -Equatorial  (SHCVS  CSBIT  INCLINED,  USER  INPUT* 

1  LONGITUDE  TO  VIEV  AT)1 

1  -GROUND  TRACK1 

I  t 

'THE  USER  IS  THEN  ASKED  TO  CHOOSE  ONE  OF  THE  FOLLOWING: ' 
'  -Unperturbed  ORBITS' 

'  -Perturbed  ORBITS' 

'  -VELOCITY  CHANGES' 

'THE  U5ER"S  CHOICE  WILL  BE  USED  IN  DEVELOPING  THE* 

|  GRAPHICAL  OUTPUT.  ' 

'THE  USER  IS  THEN  GIVEN  THE  FOLLOWING  CHOICES:  ' 

'  -CLEAR  ALL  THE  PREVIOUS  ORBITS' 

'  -CHANGE  THE  INITIAL  PARAMETERS' 

'  -CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  Anomaly' 

'  -PLOT  ANOTHER  VIEW  OF  THE  SAME  ORBIT' 


SUBROUTINE  OPTION' IOPT1 , I0PT2 , SUM , RI RAY , R JRAY , RKRAY , RARAY , 

+  TARAY,AINRAY,APRAY,TIMRAY) 

'  THIS  SUBROUTINE  GIVES  THE  USER  A  CHOICE  OF  OPERATIONS  THAT  CAN  BE 
’  PERFORMED  ON  THE  PROGRAM  AND  RETURNS  THE  USERS  CHOICE  WITH 

'  VARIABLES  IOPT1  AND  IOPT2 

DIMENSION  RIRAY(SOO) .RJRAY(SOO) ,RKRAY(500) ,RARAY(500) ,TARAY(500) , 
+  AINRAYC500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* l.YORN 
IOPT1  «  0 

'  PROMPT  USER  FOR  OPTION 

103  PRINT'S ’WHICH  OF  THE  FOLLOWING  OPTIONS  WOULD  YOU  LIKE:  ' 

PRINT*,'  1. -CALCULATE  THE  NEXT  ORBIT  USING  THE  SAME' 

PRINT'S*  PARAMETERS' 

PRINT*,'  2. -CHANGE  THE  INITIAL  PARAMETERS  OF  THE  ORBIT' 

PRINT'1','  3. -CHANGE  THE  VELOCITY  AT  A  POINT  IN  THE  ORBIT' 

PRINT'S'  (THE  UNPERTURBED  ORBIT  WILL  BE  USED)' 

PRINT'S'  4. -PLOT  ANOTHER  VIEW  OF  THE  ORBIT(S)' 

PRINT*,'  5. -QUIT' 

PRINT*, 'ENTER  i,  2,  3,  4,  OR  5:' 

READ*,I0PT2 

PRINT*,  IGPT2 

CALL  EXCMS('CLRSCRN') 

IF  (  IQPT2  .GT.  5)  THEN 
GOTO  103 
END  IF 


*  Prompt  USER  FOR  TYPE  OF  ORBIT  DESIRED 
105  IF  ( I0PT2  .EQ.  1)  THEN 

PRINT'S 'WHICH  TYPE  OF  ORBIT  WOULD  YOU  LIKE  TO  SEE,' 
PRINT'S  '  1. -Unperturbed  ORBITS' 


CRB03320 

QRB03330 

CRB03340 

0RB03350 

GRB03360 

0RB03370 

0RB03380 

0RI03390 

0RB03400 

0RB03410 

0RB03420 

0RB03430 

0X103440 

0RB03450 

0X103460 

0X103470 

0X103480 

0RB03490 

0RB03500 

CRB03510 

0RB03520 

0R103530 

0X103540 

ORB035S0 

0R103560 

ORB03570 

0RB03560 

0R903590 

0RB03600 

0RB03610 

0RB03620 

0RB03630 

0RB03640 

0RB03650 

0RB03660 

0RB03670 

ORB03680 

ORB03690 

0RB03700 

0RB03710 

0RB03720 

0RB03730 

0RB03740 

ORB03750 

ORB03760 

ORB037/0 

CRB03780 

0RB03790 

0RB038Q0 

0RBC3310 

0RB03820 

0RB03830 

0RB03840 

0RB03850 

0RB03S60 

0RB03870 


32 


PRINT*, '  2. -?ercurb«d  ORBITS' 

'PINT*, 1  ENTER  1  OR  2: ' 

REAS* , IOPTI 

PRINT*, ICPT1 

CALL  EXCMSC  'CLRSCRN' ) 

I?  KICFTi  . NE.  I)  .AND.  (IOPT1  .NE.  2;)  THEN 
PRINT1', '  INVALID  ENTRY!  ’ 

GOTO  105 
ENDIF 
END  IF 

*  PROMPT  USER  TO  CLEAR  FRSVIOUS  ORBITS 

107  IF  ( ( I0PT2  .EQ.  1)  ,0R.  (I0PT2  . EQ.  3))  THEN 

PRINT* DO  YOU  WANT  TO  CLEAR  THE  PREVIOUS  ORBITS?' 

PRINT*,' ENTER  "Y"  OR  "N"  :  ' 

READ* , YORN 

PRINT*,  YORN 

CALI,  ENCMS('Clr*crn') 

IF  ( YORN  .  EQ.  *Y')  THEN 
RIRAY(l)  *  RIRAY(NUM) 

RJRAY(l)  *  RJRAY(NUM) 

RKRAY(l)  ■  RKRAY(NUM) 

RARAY(l)  *  RARAY(NUM) 

TARAY(l)  -  TARAY(NUM) 

AINRAY(l)  -  AINRAY(NUM) 

APRAY(l)  *  APRAY(NUM) 

TIMRAY(l)  *  TIMRAY(NUM) 

SUM  -  1 

ELSEIF  (YORN  .NE.  'N')  THEN 
PRINT*, ' INVALID  ENTRY! !  ’ 

PRINT*, 'ALL  INPUTS  MUST  BE  CAPITOL  LETTERS' 

GOTO  107 

ENDIF 

ENDIF 

*  CHECK  FOR  INVALID  OPTION 

IF  ((I0PT2  .NE.  1).  AND.  (I0PT2  .NE.  2).  AND.  (I0PT2  .NE.  3)  .AND. 

+  ( I0PT2  .NE.  4). AND.  (I0PT2  .NE.  5))  THEN 

PRINT*, 'INVALID  ENTRY! ' 

GOTO  103 
ENDIF 
RETURN 
END 

*  COORDINATE  TRANSFORMATIONS 


SUBROUTINE  PQWIJKC  LAN , AP , INC , P , Q , W , I , J , K) 

THIS  SUBROUTINE  TRANSFORMS  PQV  COORDINATES  TO  IJK  COORDINATES 

DOUBLE  PRECISION  INC,P,Q,V,I,J,K,RU,R12,R13,R21,R22,R23, 

+  R31  R32  R33  LAK  AP 

Rll  *  DCOS( LAN)*DCOS( AP)  -  DSIN(LAN>*DSIN(AP)*DCOS(INC) 

R12  =  -DCOS(LAN)*DSIN(AP)  -  DS I N ( LAN )*DCOS(AP)*DCQS(INC) 

R13  =  DSIN(LAN)*DSIN( INC) 


0RB03880 

CRB03890 

DRB03900 

ORB03910 

0RB03920 

CRB03930 

ORB039&0 

0RB0395G 

0RB03960 

0RB03970 

ORB039B0 

0RB03990 

ORBOAOOO 

0RB04010 

ORBH020 

ORB04O3O 

0RB04040 

0RB04050 

0RB04060 

0RB04070 

0RB04080 

ORB04090 

0RB04100 

ORB04110 

0RB04120 

ORB04130 

ORB04140 

0RB04150 

ORB04160 

ORB04170 

0RB04180 

0RB04190 

ORB04200 

0RB04210 

0RB04220 

0RB04230 

0RB04240 

0RB04250 

0RB04260 

0RB04270 

0RB04280 

0RB04290 

0RB04300 

CRB04310 

ORB04320 

0RB04330 

0RB04340 

0RB04350 

0RB04360 

ORB04370 

ORB04380 

0RB04390 

0RB04400 

ORB04410 

ORB04420 

0RB04430 
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321  -  DS INC LANV'DCOSCAP)  +  DCOSCLAN)*DSIN(AP)*DCOS(INC) 

R22  -  -DSINC LAN>*DSINC AP)  +  DCOSC LAN)*DCOSCA?)*DCOSC INC) 

323  *  -DCOSC  LANV'DSINC  INC) 

R31  *  DSINC  A?)*DSINC  INC) 

332  *  DCCSlA?)*DSINCINC) 

333  *  DCOSC INC) 

I  »  R11*P  +  R12*Q  +  R13*V 
J  •  K21*?  +  R22*Q  +  R23*W 
X  «  R31*P  +  R32*Q  +  R33*W 
RETURN 
END 

*****************************  MM  MMMMMMMMMMMMMMMM 

SUBROUTINE  IJKPqW(LAN,AP,INC,If J,K,P,Q,W) 

♦  THIS  SUBROUTINE  TRANSFORMS  IJK  COORDINATES  TO  PQV  COORDINATES 

DOUBLE  PRECISION  INC ,  1 , J , K , P ,Q , W , R1 1 , R12»R13,R21,R22,R23, 

+  R31 ,R32,R33,LAN,A? 

Rll  »  DCOSC LAN)"DCOS(AP)  -  DSIN(LAN)*DSIN(AP)*DCOS(INC) 

R21  -  -DCOS(LAN)*DSIN(AP)  -  DSIN(LAN)*DCOS(AP)*DCOS(INC) 

R31  «  DSINC LAN)*DSIN( INC) 

R12  -  DSINC  LAN)*DCOS(AP)  +  DCOSC LAN)*DSIN(AP)*DCOSC INC) 

R22  ■  -DSINC LAN)*DSIN(AP)  +  DCOSC  LAN)*DCOSCAP)*DCOS( INC) 

R32  -  -DCOSCLAN)*DSINCINC) 

R13  ■  DSINC AP)*DSINC INC) 

R23  »  DCOSC AP)*DSINC INC) 

R33  -  DCOSCINC) 

?  *  311*1  +  R12*J  +  R13*K 
Q  «  R21*I  +  R22*J  +  R23*K 
V  »  331*1  +  R32*J  +  R33*K 
RETURN 
END 

MMM»M>>'t**i>VlM>**Vf*i>*****M  ************************************ 

SUBROUTINE  IJKRSWC LAN , AL , INC , I , J , K , R , S , V) 

*  THIS  SUBROUTINE  CHANGES  FROM  IJK  COORDINATES  TO  RSV  COORDINATES 

DOUBLE  PRECISION  INC,I>J,K#R>S,W,Rll,R12>R13,R21,R22,R23f 
+  R31,R32,R33,LAN,AL 

Rll  -  DCOSC LAN)*DCOS( AL)  -  DSINC LAN)*DCOS( INC)*DSIN(AL) 

R12  *  DSINC LAN)*DCOSCAL)  +  DSIN(AL)*DCOSCLAN)*DCOSCINC) 

R13  ■  DSINC INC)*DSIN(AL) 

R21  »  -DCOSC LAN)*DSINCAL) -DSINC LAN)*DCOS(INC)*DCOSCAL) 

R22  ■  -DSINC LAN)*BSINCAL)  +  DCOSC  LAN)*DCOSC INC)*DCOSC AL) 

R23  ■  DSINC INC)*DCOSCAL) 

R31  »  DSINC LAN)*DSIN( INC) 

R32  «  -DCOSC  LAN)*DSINC INC) 

R33  *  DCOSCINC) 

R  -  R11*I  +  R12*J  +  R13*K 
S  «  R21*I  +  R22*J  +  R23*K 
W  *  R31*I  +  R32*J  +  R33*K 
RETURN 
END 


ORB04A40 

0RB044S0 

ORBOAAAO 

0RB0AA70 

0RB0AA60 

0RB0AA90 

0RB0A500 

0RB0A510 

ORBOAS20 

0RB0A530 

ORBOASAO 

ORBOASSO 

0RB04560 

0RB0A570 

ORBOASSO 

0RB0AS90 

ORBOA6O0 

ORB04610 

0RB0A620 

0RB0A630 

ORBOASAO 

0RBQA650 

ORB0A66O 

ORB0A67O 

0RB0A680 

0RB0A690 

ORB0A700 

0RB0A710 

0RB0A720 

0RB0A730 

ORBOA7AO 

ORBOA7SO 

ORBOA760 

0RB0A770 

ORBOA7SO 

0RB0A790 

ORB0A800 

ORB0A810 

0RB0A820 

ORBOA830 

ORBOASAO 

0RB0A850 

ORBOA860 

0RB0A870 

0RB0A880 

ORB0A890 

0RB0A900 

0RB0A910 

0RB04920 

0RB0A930 

ORBOA9AO 

ORB0495O 

0RB04960 

ORB04970 

OR304980 
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SUBROUTINE  RSV  X JK( LAN , AL, INC , R , S , V , I , J , K) 

THIS  SUBROUTINE  CHANGES  FROM  RSV  COORDINATES  TO  IJK  COORDINATES 

DOUBLE  PRECISION  INC,R,S,W,I,J.K,RlllR12fR13,R21,R22,R23, 

+  R31,R32,R33,LAN,AL 

Rll  «  DCOS(LAN)*DCOS(AL)  -  DSIN(LAN)*DCOS(INC)*DSIN(AL) 

R2I  ■  DSXNC LAN)*DCOS(AL)  +  DSIN(AL)*DCOS(LAN)*DCOS(INC) 

R31  ■  DSIN( INC)*DSIN(AL) 

R12  «  -DCOSC LAN)*DSIN( AL) -DSIN( LAN)*DCOS( INC)*DCOS( AL) 

R22  *  -DSIN(LAN)*DSIN(AL)  +  DCQS(LAN)*DCOS(INC)*DCOS(AL) 

R32  »  DSIN( INC)*DCOS(AL) 

R13  «  DSIN(LAN)*DSIN(INC) 

R23  *  -DCOSC LAN)*DSIN( INC) 

R33  «  DCOS( INC) 

I  «  RII*R  +  R12*S  +  R13*W 
J  «  R21*R  +  R22*S  +  R23*W 
K  «  R31*R  +  R32*S  +  R33*V 
RETURN 
END 


SUBROUTINE  PQWRSW(TA , P , Q , V , R , S , VN) 

THIS  SUBROUTINE  CHANGES  FROM  PQW  COORDINATES  TO  RSV  COORDINATES 

DOUBLE  PRECISION  P,Q,W,R,S,VN,RllfR12,R13,R21lR22,R23, 

+  R31,R32,R33,TA 

Rll  *  DCOSC TA) 

R12  -  DSIN(TA) 

R13  *  0.  0 
R21  «  -DSIN(TA) 

R22  *  DCOSC TA) 

R23  -  0.  0 
R31  *  0.0 
R32  -0.0 
R33  *  1.0 

R  ■  RU*P  +  R12*Q  +  R13*V 
S  =  R21*P  +  R22*Q  +  R23*V 
VN  «  R31*P  +  R32*Q  +R33*V 
RETURN 
END 


SUBROUTINE  RSWPQV( TA , R , S , V , P , Q , WN) 

*  THIS  SUBROUTINE  CHANGES  FROM  RSV  COORDINATES  TO  PQW  COORDINATES 

DOUBLE  PRECISION  R,S,V,P,Q,WN,R11,R12,R13,R21,R22>R23, 

+  R31,R32,R33,TA 

Rll  *  DCOS( TA) 

R21  *  DSINCTA) 

R31  =  0.0 
R12  =  -DSINCTA) 


ORB04990 

DRS05000 

0RB05010 

0RB05020 

0RB05030 

0RB05040 

0RB05050 

ORBO506O 

0RB05070 

ORB05080 

ORB05090 

ORBOSIOO 

0RI05110 

ORB05120 

0RB05130 

ORB05140 

0RB0S150 

0RB03160 

ORB0517O 

0RB05180 

0RB05190 

0RB05200 

0RB05210 

ORB05220 

ORB05230 

0RB03240 

0RB05250 

ORB05260 

ORB05270 

ORB05280 

ORB05290 

OP.B05300 

ORB05310 

ORB05320 

ORB05330 

GRB05340 

ORB05350 

ORB05360 

ORB05370 

ORB05380 

0RB05390 

ORB05400 

0RB05410 

ORB05C20 

ORB05430 

ORB05440 

ORB05450 

ORB05460 

ORB05470 

ORB05480 

ORB05490 

ORB05500 

0RB05510 

ORB05520 

ORB05530 

ORB05540 


R22  «  DCOS(TA) 

R32  *  0. 0 
R13  *  0.0 
R-3  «  0.0 
R33  *  1.0 

P  *  R11*R  +  R12*S  +  R13*V 
Q  «  R21*R  +  R22*S  +  R23*V 
VN  «  R31*R  +  R32*S  +R33*W 
RETURN 
END 


STORE  ELEMENTS  IN  ARRAYS 


SUBROUTINE  STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY,TARAY,NUM, 
+  I,AP,AINRAY,APRAY,TT,TIMRAY) 

THIS  SUBROUTINE  STORES  THE  POSITION  AND  ELEMENTS  IN  ARRAYS  IN 
SINGLE  PRECISION  FORM,  FOR  PLOTTING 

DOUBLE  PRECISION  RI ,RJ,RK,R,TA,I ,AP,TT 

DIMENSION  RIRAY(SOO) , RJRAY(SOO) ,RKRAY(500) ,RARAY(500) ,TARAY(500) , 
+  AINRAY(500),APRAY(500),TIMRAY(500) 

RIRAYCNUM)  «  SNGL(RI) 

RJRAY(NUM)  *  SNGL(RJ) 

RKRAY(NUM)  ■  SNGL(RK) 

RARAY(NUM)  *  SNGL(R) 

TARAY(NUM)  -  SNGL(TA) 

AINRAY(NUM)  ■  SNGL(I) 

APRAY(NUM)  -  SNGL(AP) 

TINRAY(NUM)  «  SNGL(TT) 

RETURN 

END 


*  INITIAL 


POSITION,  VELOCITY 


SUBROUTINE  INPUTS(RI,RJ,RK,R,VI,VJ,VK,V,MU,QUIT  FT) 

THIS  SUBROUTINE  GIVES  THE  USER  A  CHOICE  TO  EJTTP  .'ITER  THE 
'  INITIAL  POSITION  AND  VELOCITY  VECTOR  OR  TO  l'£  ftCi  PROGRAM 
'  CALCULATE  THE  INITIAL  POSITION  AND  VELOCITY  If  R  PROMPTED 

f  INPUTS 

'  SUBROUTINES  CALLED  FROM  THIS  SUBROUTINE: 

'  INELTS  *  Prompts  USER  FOR  ORBITAL  ELEMENTS 

f  IPOS  *  PROMPTS  USER  FOR  INITIAL  POSITION  (IJK) 

IVEL  »  PROMPTS  USER  FOR  INITIAL  Velocity  (IJK) 

DOUBLE  PRECISION  RI,RJ,RK,R,VI,VJ,VK,V,MU,PI 
CHARACTER*^  QUIT 

f  PROMPT  USER  FOR  METHOD  TO  ENTER  INPUTS 

195  PRINT*, 'IN  WHICH  MANNER  WOULD  YOU  LIKE  TO  INPUT  THE  INITIAL' 


0RB05550 

0RB05560 

0RB03570 

0RB05560 

0RB05590 

0RB05600 

0RB0S610 

0RB05620 

0RB05630 

ORB05640 

0RB05650 

0RB05660 

ORB05670 

0RB05680 

ORB05690 

0RB05700 

OPB05710 

0RB05720 

0RB05710 

0RB05740 

0RB05750 

0RB05760 

0RB0S770 

ORB05780 

0RB0S790 

0RB0S800 

ORB05810 

0RB05820 

ORB05830 

ORB05840 

0R305850 

ORB05860 

0RB05870 

ORB05880 

0RB0589C 

ORB05900 

0RB05910 

ORB05920 

ORB05930 

ORB05940 

0RB05950 

ORB05960 

0RB05970 

ORB05980 

ORB05990 

ORB06000 

0RB06010 

0RB06020 

0RB06030 

ORB06040 

0RB06050 

0RB06060 

ORB06070 

ORB06080 

ORB06090 

0RB06100 
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PRINT*,' POSITION  AND  VELOCITY  OF  THE  SATELLITE?' 

PRINT*, '  I:  BY  Inputting  THE  INITIAL  POSITION  AND  VELOCITY' 
PRINT*,'  VECTORS  IN  THE  PERIFOCAL  COORDINATE  SYSTEM  (IJK)' 

PRINT'','  2:  5Y  LETTING  THE  SATELLITE  BE  PLACED  ON  THE  "I"' 
PRINT*,’  AXIS  OF  THE  (IJK)  SYSTEM  AT  A  DESIRED  RADIUS  OF' 

PRINT'S'  PERIGEE(RP)  A»..  INPUTTING  EITHER  A  DESIRED  RADIUS 

PRINT'S'  OF  APOGEE(RA) ,  A  DESIRED  ECCENTRICITY(E) ,  OR  THE' 

PRINT*, '  DESIRED  VELOCITY  AT  THAT  RADIUS,  AND  A  DESIRED' 

PRINT* , '  INCLINATION^  I) .  ' 

PRINT*,’  3:  QUIT' 

PRINT'S 'ENTER  1,  2  OR  3: ' 

READ*, ICHC 

PRINT*, ICHC 

CALL  EXCMS('CLRSCRN') 

*  USER  INPUTS  POSITION  AND  VELOCITY  VECTORS 
IF  (ICHC  .EQ.  1)  THEN 

CALL  ir«JS(RI ,RJ,RK,R) 

CALL  I VEL( VI , V J , VK , V , R , MU ) 

*  USER  INPUTS  ORBITAL  ELEMENTS  TO  GET  POSITION  AND  VELOCITY 
ELSEIF  (ICHC  .EQ.  2)  THEN 

CALL  INELTS(RI,RJ,RK,R,VI,VJ,VK,V,MU,PI) 

*  STOP  PROGRAM 

ELSEIF  (ICHC  .EQ.  3)  THEN 
QUIT  »  'X' 

ELSE 

PRINT’S ' INVALID  ENTRY!  TRY  AGAIN! ' 

GOTO  195 
ENDIF 
RETURN 
END 


iV>Y«iVVn>iV.VVnVVn>iVVnVfti>iVVrjV*i>iV)ViViVAvV#iV^VWc#*iVi>**iYiV^*rt**rttoV*>>*VtfoMr*ViYilnflr»lr 


SUBROUTINE  IFOS(RI ,RJ,RK,R) 

*  THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  INITIAL  POSITION  OF  THE 

*  Satellite  IN  GEOCENTRIC-EQUATORIAL  COORDINATE  SYSTEM 

DOUBLE  PRECISION  RI,RJ,RK,R 

CHARACTER* 1,  CHOICE 
LOGICAL  CORREC 
COP.REC  »  .FALSE. 

*  PROMPT  USER  FOR  VELOCITY  VECTOR 
180  IF(.  NOT.  CORREC)  THEN 

CALL  EXCMS('CLRSCRN') 

PRINT'S 'ENTER  RADIUS  VECTOR  VALUES  IN  "KM"' 

PRINT*, 'RADIUS  OF  THE  EARTH  =  6400  KM' 

CORREC  =  .  TRUE. 

PRINT’’', 'ENTER  RI  : ' 

RE AD*, R I 

PRINT’S 'RI  =  P.1,1  KM' 

PRINT*, 'ENTER  RJ  : ' 


ORB06110 

0R306120 

ORB06130 

ORB06140 

0RB06150 

ORB06160 

0RB06170 

0RB06160 

0RB06190 

0RB06200 

0RB06210 

0RB06220 

0RB06230 

0RB06240 

0RB06250 

ORB06260 

0RB06270 

GRB06280 

0RB06290 

0RB06300 

0RB06310 

0RB06320 

0RB06330 

0RB06340 

0RB06350 

0RB06360 

ORB06370 

0RB06380 

0RB06350 

ORB06400 

0RB06410 

0RB06420 

0RB06430 

0RB06440 

ORB06450 

0RB06460 

0RB06470 

0RB06480 

0RB06490 

ORB06500 

0RB06510 

ORB06520 

0RB06530 

ORB06540 

ORB06550 

0RB06560 

0RB06570 

0RB06580 

0RB06590 

0RB06600 

0RB06610 

0RB06620 

0RB06630 

ORB06640 

0RB06650 

0RB06660 


HEAD'S  RJ 

PRINT* , 1 RJ  «  1  .RJ.'KM* 
f HINT'S* ENTER  RK  :  1 
K&nL<M ;  R.\ 

FAINT'S*  HK  =  1  ,RK,  *KM' 

CALCULATE  TOTAL  R 

R  ■  DSQRT((RI**2)  +  (RJ**2)  +  (RK**2)) 

PRINT* , 1 R  »  ’ ,R, ‘KM* 

IF  (R  .LE.  6400.0)  THEN 

PRINT’S  'RADIUS  TO  SMALL! !  ESTER  NEW  VALUES! !  ' 
GOTO  ISO 

ENDIF 


CHECK  WITH  USER  THAT  Valuas  ARE  CORRECT 
PRINT*, 'ARE  THESE  VALUES  CORRECT?' 
PRINT’S'  ENTER  "Y"  OR  "N" 

READ*, CHOICE 
CHOICE  «  ' Y' 

PRINT’S  CHOICE 
IF  (CHOICE. EQ. ' Y' )  THEN 
CORREC  «  .TRUE. 

ENDIF 
GOTO  180 
ENDIF 
RETURN 
END 


i’fVh>*i'n'n'n,nV,ViViV*i'n'n’M,nV*,VA***»V**Vr*AWnV*in>*#ir*i'n¥^*»>*iVA**VhVAV,nV*iWM'nWn1nViKf 

SUBROUTINE  IVEL( VI , Vi , VK , V , R ,MU) 

*  THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  INITIAL  VELOCITY'  OF  THE 

*  S’cellice 

DOUBLE  PRECISION  VI ,VJ,VK,V,R,VCIR,VMAX,MU 

CHARACTER* 1 ,  CHOICE 
LOGICAL  CORREC 
CORREC  «  .FALSE. 

*  CALCULATE  ESCAPE  VELOCITY  AND  CIRCULAR  VELOCITY  AND  PROMPT  USER 

*  FOR  VELOCITY  VECTOR 
190  IF(. NOT. CORREC)  THEN 

CALL  EXCMS( 'CLRSCRN' ) 

VCIR  *  DSQRT(MU/R) 

VMAX  »  DSQRT((2. 0*MU)/R) 

PRINTS 'CIRCULAR  VELOCITY  «  ', VCIR  'KM/SEC' 

PRINT'S 'MAXIMUM  VELOCITY’  «  ' , VMAX,  KM/SEC' 

CORREC  *  . TRUE. 

PRINTS 'ENTER  VELOCITY  VECTOR  IN  (KM/SEC)' 

PRINT*, 'ENTER  VI  : ' 

READ*, VI 

PRINT'S 'VI  =  ', VI, 'KM/SEC' 

PRINTS' ENTER  VJ  :  ' 

READ* , VJ 


ORB06670 

0RB06680 

0RB06690 

0RB06700 

0RB06710 

0RB06720 

0RB06730 

0RB06740 

ORR06750 

ORB06760 

0RB06770 

0RB06780 

0RB06790 

ORB06800 

0RB06810 

ORB06820 

ORB06830 

ORB06840 

ORB06850 

0RB06860 

0RB06870 

0RB06880 

ORB06890 

0RB06900 

0RDQ6910 

0RB06920 

ORB06930 

ORB06940 

ORB06950 

OR306960 

0RB06970 

0RB06980 

ORB06990 

0RB07000 

ORB07010 

ORB07020 

ORB07030 

0RB07040 

ORBO7O50 

ORB07060 

0RB07070 

ORB07080 

ORB07090 

ORB07100 

ORB07110 

0RB07120 

ORB07130 

ORB07140 

ORB07150 

ORB07160 

ORB07170 

0RB07180 

ORBO7190 

ORB07200 

ORB07210 

0R307220 


PRINT'S 'VJ  *  *,VJ, 'KM/SEC' 

PRINT'S 'ENTER  VK  : ' 

READ*,VK 

PR  INT* , '  VK  -  ' ,  VK , '  KM/SEC ' 

CALCULATE  TOTAL  VELOCITY  (V) 

V  «  DSQRT((VI**2)  +  (VJ**2)  +  (VK**2)) 

PRINT* , 1 V  »  ‘,V,f KM/SEC' 

CHECK  WITH  USER  THAT  VALUES  ARE  CORRECTS 
PRINT*, 'ARE  THESE  VALUES  CORRECT?' 

PRINT*, '  ENTER  "Y"  OR  "N" 

READ*, CHOICE 
CHOICE  «  'Y' 

PR INT*, CHOICE 
I?  (CHOICE. EQ. 'Y')  THEN 
CORREC  «  .TRUE. 

END  IF 

IF  (V  .  5E.  VMAX)  THEN 

PRINT*, 'VELOCITY  IS  GREATER  THAN  THE  ESCAPE  VELOCITY!!’ 
PRINT*, 'RE-ENTER  VELOCITY! I ! ' 

CORREC  *  .FALSE. 

ENDIF 
GOTO  190 


ENDIF 

RETURN 

END 


SUBROUTINE  INELTS(RI.  RJ,  RK,  R,  VI,  VJ,  VK,  V,MU,PI) 
SATELLITE  PLACED  ON  fl'  AXIS  AND  USER  SUPPLY  ORBITAL  ELEMEJ 
GET  INITIAL  POSITION  AND  VELOCITY 


ELEMENTS  TO 


DOUBLE  PRECISION  RI,RJ,RK,R,VI,VJ,VK,V,MU,I,ENR,A,E,RP,RA,PI,VMAX 
CHARACTER* 1 .CHOICE 

PROMPT  USER  FOR  PERIGEE  RADIUS 

PRINT*, 'ENTER  RADIUS  OF  PERIGEE(RP)  IN  (KM),  FOR  EXAMPLE:' 

PRINT*, 'LOW  EARTH  ORBIT  (LEO),  RP  ■  6600.0  KM' 

PRINT*, 'GEOSYNCROUNOUS  ORBIT,  RP  *  42241.1  KM' 

PRINT*, 'ENTER  RP:  ' 

PRINT*, '"RP"  MUST  BE  >  6400KM' 

READ*,RP 
PRINT*, RP 

CHECK  FOR  VALID  RADIUS 
IF  (RP  .LT.  6400.0)  THEN 

PRINT*, 'YOUR  "RP"  IS  TO  SMALL!!' 

GOTO  198 
ENDIF 

PROMPT  USER  FOR  TYPE  OF  INPUT 

PRINT'S 'DO  YOU  WANT  TO  ENTER  THE  ECCENTRICITY  (E),  ' 

PRINT'S ' RADIUS  OF  APOGEE  (RA).  OR  VELOCITY  (V)?’ 

PRINT'S 'ENTER  "E",  "R",  OR  "V":  ' 


0RB07230 

0RB07240 

0RB07250 

ORB07260 

ORB0727O 

0RB07280 

ORB07290 

0RB07300 

0RB0731O 

ORB07320 

0RB07330 

0RB07340 

0RB073S0 

ORB07360 

ORB07370 

0RB07380 

0RB07390 

ORB07400 

0RB07410 

0RB07420 

0RB07430 

0RB07440 

ORB074S0 

0RB07460 

0RB07470 

ORB07480 

ORB07490 

ORB07500 

0RB07510 
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READS  CHOICE 
PRISTS  CHOICE 
CALL  EXC!!S(  'CLRSCRN  ) 

*  USER  ENTERS  Eccentricity  AND  SEMI-MAJOR  AXIS,  ENERGY  AND  VELOCITY 

*  IS  CALCULATED  IN  THAT  ORDER 
IF  (CHOICE  .EQ.  'E')  THEN 

PRINT*, 'ENTER  ECCENTRICITY  (E):  ’ 

PRINT*,' 0.0  <■  E  <  1.0' 

READ*,E 
PRINT'S  E 

*  CHECK  FOR  VALID  ECCENTRICITY 

IF  <(E  .LT.  0.0)  .OR.  (E  .GE.  1.0))  THEN 
PRINT*, 'INVALID  "E"' 

GOTO  198 
ENDIF 

A  ■  RP/U-E) 

ENR  -  -MU/( 2. 0*A) 

V  «  DSQRT(2*(ENR+(MU/RP) ) ) 

*  USER  INPUTS  RADIUS  OF  APMEE  AND  ECCENTRICITY  IS  CALCULATED 

*  THEN  SEMI -MAJOR  AXIS,  ENERGY  AND  THEN  VELOCITY. 

ELSEIF  (CHOICE  .  EQ.  'R')  THEN 

PRINT’S 'ENTER  RADIUS  OF  APOGEE  (RA)  IN  KM:' 

PRINTS' "RA"  MUST  EE  >«"RP",  "RP"  ■  '  ,R? 

READ*, RA 
PRINT*, RA 

*  CHECK  FOR  VALID  RADIUS  OF  APOGEE 
IF  (RA  .LT.  RP)  THEN 

PRINT'S  'YOUR  "RA"  IS  TO  SMALL!!' 

GOTO  198 
ENDIF 

E  «  (RA-RP)/(RA+RP) 

A  *  P.P/O-E) 

ENR  «  -MU/(2. 0*A) 

V  ■  DSQRT(2*(ENR+(MU/RP))) 

*  USER  INPUTS  MAGNITUDE  OF  VELOCITY,  PROGRAM  PROVIDES  CIRCULAR 

*  AND  ESCAPE  VELOCITY  FOR  COMPARISON  AND  TO  CHECK  FOR  VALID 

*  INPUTS 

ELSEIF  (CHOICE  .EQ.  ' V* )  THEN 

PRINT*, 'ENTER  VELOCITY  IN  KM/SEC: ' 

PRINT*, 'THE  MINIMUM  VELOCITY  ALLOWED  IS  FOR  A  CIRCULAR  ORBIT' 
VC IRC  -  SQRT(SNGL(MU/RP)) 

PRINT*, 'ORBIT.  V(Circular)  ■  \VCIRC,'  KM/S' 

VMAX  -  DSQRT( 2*( MU/RP) ) 

PRINT*, 'THE  MAXIMUM  VELOCITY  <  ',VMAX,'  KM/S' 

READ*,V 
PRINT*, V 

IF  (V  .LT.  VCIRC)  THEN 

PRINT*, 'VELOCITY  TO  SMALL!' 

GOTO  198 

ENDIF 

IF  (V  .GE.  VMAX)  THEN 
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PRINT*,' VELOCITY  70  GREAT! ! 1 
GOTO  1*38 
END  IF 

PRINT*, ' INVALID  ENTRY!  TRY  AGAIN* 

GOTO  198 
END  IF 

INCLINATION  NEEDED  TO  GIVE  Velocity  A  Direction 
PRINT’S 'ENTER  INCLINATION  (I)  IN  DEGREES:' 
READ*, I 
PRINT’S  I 

I  «  (PI/180. 0)*I 
VX  -  V*DSIN(I) 

VJ  *  V*DCOS( I) 

VI  «  0.0 

RADIUS  VECTOR  SET 

RI  «  RP 

RJ  »  0.  0 

RX  «  0.  0 

R  *  RP 

RETURN 

END 


*  CALCULATE  THE  ORBITAL  ELEMENTS 

SUBROUTINE  CALCELC RI , R J , RK , R , VI , VJ , VX , V , El , E J , EK , E , A , I , LAN , 

+  LPtTA,PER,EA,MA,AP,AL,TF,P,PI ,MU,MM,N,H,HI ,HJ) 

*  THIS  SUBROUTINE  CALLS  THE  INDIVIDUAL  SUBROUTINES  TO  CALCULATE  THE 

*  ORBITAL  ELEMENTS 

*  THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES( RETURNED  VALUES) 


ENERGY 

ANGMOM 

NODE 

LATREC 

ECC 

SMAXIS 

INCL 

ASNODE 

ARP 

IJKPQW 

TANOM 

ARLAT 

LONPER 

TLON 

PERIOD 

ECCAN 

MEANMO 

HEANAK 

TFLGHT 


ENERGY  PER  MASS  (ENR) 

ANGULAR  MOMENTUM  (H,HI,HJ,HK) 

NODE  VECTOR  (N,NX,NJ,NK) 
SEMI-LATUS  RECTUS  (P) 

ECCENTRICITY  (E,EI,EJ,EK) 
SEMI-MAJOR  AXIS  (A) 

INCLINATION  (I) 

LONGITUDE  OF  ASCENDING  NODE  (LAN) 
ARGUMENT  OF  PERIGEE  (AP) 

' IJK*  SYSTEM  TO  ' PQW '  SYSTEM 
TRUE  ANOMALY  (TA) 

ARGUMENT  OF  LATITUDE  (AL) 
LONGITUDE  OF  Perigee  (LP) 

TRUE  LONGITUDE  (TL) 

PERIOD  (PER) 

ECCENTRIC  ANOMALY  (EA) 

MEAN  MOTION  (MM) 

MEAN  ANOMALY  (MA) 

TIME  OF  FLIGHT  (TF) 


DOUBLE  PRECISION  RI,RJ,RK,R,VI ,VJ,VK,V,EI ,EJ,EX,E,A,I ,LAN,AL, 
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■f  IT  ,TA ,  PER .  SA ,  MA ,  A? ,  TF ,  HI ,  H  J ,  HjC ,  H ,  X I ,  N*J ,  NX .  X ,  P ,  ?I ,  MU ,  MM ,  ENR , 

+  TL,S?,H3,RV,NP,NQ,NW 


CALL  ENERGYC  V ,  R  ,MU ,  ENR) 

CALL  ASGMQiK RI , R J , RK , VI , VJ , VK ,HI, HJ , KK , H) 

CALL  N:CE(HI ,HJ,NI ,NJ,NK,N) 

CALL  LATREC(H,?,HU) 

CALL  ECC(RI,RJ(RK>R,VI»VJ,VK,V,EIIEJ,EK,E,MU) 

CALL  SMAXIS(MU,ENR,A) 

CALL  INCL(HK,H,I,PI) 

*  SPECIAL  CASE  IF  INCLINATION*  «  0.0 
IF  (I.NE.O.O)  THEN 

CALL  ASSQDE( NI ,  N , LAN , NJ , PI ) 

CALL  ARP(NI ,NJ,N,EI,EJ,EK,E,AP,PI ,NP,NQ,LAN) 

ELSE 

LAN  -  0. 0 
A?  *  0.  0 
SNDIF 

*  COORDINATE  TRANSFORMATION  OF  'R'  AND  'V'  VECTORS 
CALL  I JKPQV( LAN , AP , I , R I , R J , RK , RP , RQ , RV ) 

CALL  I JK PQW( LAN , AP , I , N I , N J , NK , NP , NQ ,  NV ) 

CALL  TANOMCEI.EJ.EK.E.RI.RJ.RK.RP.RQ.RV.R.VI.VJ.VK.TA.PI) 

*  SPECIAL  CASE  FOR  Inclination  «  0.0 
IF  (I  .NE.  0.0)  THEN- 

CALL  ARLAT( NI , NJ , NK , N , RI , R J , RK , S , AL , PI ,TA , AP) 

ELSE 

AL  »  TA 
END  IF 

CALL  LONPEKC  LAN , A? , LP) 

CALL  TLQN(LAN,AP,TA,TL) 

CALL  PSRIOD( A , PER , PI , MU) 

CALL  ECCAN( E ,TA , EA , PI ) 

CALL  MEANMO(A,MM,MU) 

CALL  MEANAN(EA,E,MA) 

CALL  TFLGHT( MM , MA , TF ) 

RETURN 

END 

SUBROUTINE  ENERGY( V , R ,MU , ENR) 

*  THIS  SUBROUTINE  CALCULATES  THE  ENERGY  OF  THE  ORBIT 

DOUBLE  PRECISION  V,R,MU,ENR 


ENR  =  ((V**2)/2)  -  (MU/R) 

RETURN 

END 


>'nViV*iVi'tjKc>,o*nKc*Ai'fVfiV>ViViV»lriyi>jViVVir#iHr^nWt*i'r;ViViViVVf#)VjViViVVriV***iV»WnVVr*^rA*iV*iV*iyilriV 


SUBROUTI NE  ANGMOMC R I , R J , RK , VI , V J , VK , HI , H J , KK , H ) 
THIS  SUBROUTINE  CALCULATES  THE  ANGULAR  MOMENTUM 
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DOUBLE  PRECISION  HI.RJ.RK.VI.VJ.VK.HI.HJ.HK.H 

HI  »  (RJ  *  VS)  -  (RK  *  VJ) 

HJ  *  CSK  *  vn  -  (F.I  *  VS) 

HK  «  <RI  *  VJ)  -  (RJ  *  VI) 

H  ■  DSQRT((HI**2)  +  (KJ**2)  +  (HK**2)) 

RETURN 

END 


SUBROUTINE  N'ODE(HI ,HJ,NI ,NJ,NK,N) 

THIS  SUBROUTINE  CALCULATES  THE  NODE  VECTOR 

DOUBLE  PRECISION  HI ,HJ,NI .NJ^NK.N 

NI  »  -HJ 
NJ  =  HI 
NK  «  0.  0 

N  »  DSQRT((NI**2)  +  (NJ**2)) 

RETURN 

END 


SUBROUTINE  LATREC(H,P,MU) 

*  TillS  SUBROUTINE  CALCULATES  THE  SEMI-LATUS  RECTUM 

DOUBLE  PRECISION  H,P,MU 

P  *  (H**2)/MU 

RETURN- 

END 

Vn'n>VWra'nVVn':uVn>>V*5,nVi'ri'(i'rVnWM>*Vr**i>Vf*iWfiV*iWr*Vr#i,n,rynVTV*i'f*yn,rff***A- 


SUBROUTINE  ECC(RI ,RJ,RK,R,VI ,VJ,VK,V,EI ,EJ,EK,E,MU) 

THIS  SUBROUTINE  CALCULATES  THE  ECCENTRICITY 

DOUBLE  PRECISION  RI.RJ.RK.R.V^VJ.VK.VjEI.EJ.EK.E.MU.DOT 

CALCULATE  DOT  PRODUCT  OF  'R'  AND  * V'  VECTORS 
DOT  -  (RI*VI)  +  (RJ*VJ)  +  (RK*VK) 

El  »  ( 1. OD+OO/MU)  *  ( ( ( V**2 )  -  (MU/R))  *  RI  -  (DOT)*VI) 

EJ  *  (1  OD+OO/MU)  *  (((V**2)  -  (MU/R))  *  RJ  -  (DOT)*VJ) 

EK  »  (1. OD+OO/MU)  *  (((V**2)  -  (MU/R))  *  RK  -  (DOT)*VK) 

E  «  DSQRT((EI**2)  +  (EJ**2)  +  (EK**2)) 

RETURN 

END 


SUBROUTINE  SMAXIS(MU>ENRfA) 

THIS  SUBROUTINE  Calculates  THE  SEMI-MAJOR  AXIS 
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DOUBLE  PRECIS  ICS  MU,ENR,A 

a  =  -::?;/(  2*enr) 

RETURN 

END 


SUBROUTINE  INCL(KK,H,I,PI) 

*  THIS  SUBROUTINE  CALCULATES  THE  INCLINATION 

*  V  ALWAYS  LESS  THAN  180  DEGREES 

DOUBLE  PRECISION  HK,H,I,PI 

I  «  DACOS(HK/H) 

RETURN 

END 


SUBROUTINE  ASNODE( NI , N , LAN , N J , PI ) 

*  THIS  SUBROUTINE  CALCULATES  THE  LONGITUDE  OF  THE  ASCENDING  NODE 

*  IF  ’NJ*  >  0  THEN  'LAN'  <  180  DEGREES 

DOUBLE  PRECISION  NT  ,N,LAN,NJ,?I 

LAN  *  DATAN2(NJ,NI) 

IF  (LAN  .LT.  0.0)  THEN 
LAN  «  (2*PI)  +  LAN 
END  IF 
RETURN- 
END 


SUBROUTINE  ARP( NI , NJ , N , El , E J , EK , E , AP , PI ,NP , NQ , LAN) 

*  THIS  SUBROUTINE  CALCULATES  THE  ARGUMENT  OF  Parigaa 

*  IF  'EK'  GREATER  THAN  0  THEN  'AP'  <  180 

*  VARIABLE  TEMP  USED  AS  A  Temporary  VALUE  FOR  ARCTAN 

DOUBLE  PRECISION  NT , NJ ,N , El , E J , EK ,E , AP , PI ,NQ , NP ,TEMP , LAN 

IF  C  (El  .EQ.  0.0)  .AND.  (EJ  .  EQ.  0.0))  THEN 
AP  »  0.  0 
ELSE 

TEMP  «  DATAN2(EJ,EI) 

IF  (TEMP  .GT.  LAN)  THEN 
AP  «  TEMP  -  LAN 

ELSE 

AP  =  (2*PI)  -  (LAN  -  TEMP) 

ENDIF 

IF  (  AP  .LT.  0.0)  THEN 
A?  -  (2*PI)  +  AP 
ENDIF 

IF  (AP  .GT.  (2*PI))  THEN 
A?  =  A?  -  (2"PI) 
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ENDIF 

ESDI? 

RETURN 

END 

SUBROUTINE  TANOM( El , E J , EK , E , RI . R J , RK , RP , RQ , RV , R , VI , VJ , VK , 

+  7A,?I) 

*  THIS  SUBROUTINE  CALCULATES  THE  TRUE  Anomaly 

*  IF  (R  DOT  V)  >  0  THEN  TA  <  180  DEGREES 

DOUBLE  PRECISION  DOT.EI.EJ.EK.E.R^RJ.RK.R.VI.VJ.VK.TA.PI, 
+  RP.SQ.RV 

TA  -  DATAN2(RQfRP) 

IF  (TA  .LT.O.Q  )  THEN 
TA  -  (2  *  PI)  +  TA 
ENDIF 
RETURN 
END 


SUBROUTINE  ARLAT( NT , NJ , NK , N , RI , R J , RK , R , AL , PI ,TA  t AP) 

*  THIS  SUBROUTINE  CALCULATES  THE  ARGUMENT  OF  LATITUDE 

*  IF  (RK  >  0)  THEN  AL  <  180  DEGREES 

DOUBLE  PRECISION  NT , NJ , NK , N , RI , R J , RK , R , AL, PI ,TA , A? 

AL  *  TA  +  A? 

RETURN- 

END 

SUBROUTINE  LONPER(LAN,AP,LP) 

*  THIS  SUBROUTINE  CALCULATES  7HE  LONGITUDE  OF  PERIGEE 


DOUBLE  PRECISION  LAN,AP,LP 

LP  «  LAN  +  A? 

RETURN 

END 


SUBROUTINE  TLON( LAN , AP , TA , TL) 

THIS  SUBROUTINE  CALCULATES  THE  TRUE  LONGITUDE  AT  EPOCH 

DOUBLE  PRECISION  LAN,AP,TA,TL 

TL  *  AP  +  LAN  +  TA 

RETURN- 

END 


QRB10590 

ORB 10600 

CR310610 

ORB 10620 

ORB 10630 

ORB 10640 

0RB10659 

ORB 10660 

ORB 10670 

ORB 10680 

ORB 10690 

ORBIO/OO 

0RB10710 

0RB10720 

ORB10730 

0RB10740 

ORB10750 

0RB10760 

0RB10770 

ORB10780 

0RB10790 

ORB 10800 

0RB10810 

CRB10820 

ORB  10830 

ORB 10840 

0RB10850 

0RB10860 

ORB 10870 

ORB10S80 

ORB10S90 

ORB10900 

0RB10910 

0RB10920 

ORB10930 

ORB 10940 

0RB109S0 

ORB 10960 

0RB10970 

ORB10980 

ORB10990 

ORB 11000 

ORBllOlO 

0RB11020 

0RB11030 

QRB11040 

ORBllOSO 

0RB11060 

ORB11070 

ORB 11080 

0RB11090 

0RB11100 

ORBllllO 

0RB11120 

0RB11130 


45 


AittVi 


* 


SUBROUTINE  ?ERIOD(A#PER,PI ,MU) 

THIS  SUBROUTINE  CALCULATES  THE  PERIOD 


DOUBLE  PRECISION*  A,?ER,PI,MU 

PER  «  2.  QD+0Q*(PI)*DSQRT((A**3)/MU) 
RETURN* 

END 


SUBROUTINE  ECCAN(E,TA,EA,PI) 

*  THIS  SUBROUTINE  CALCULATES  THE  ECCENTRIC  Anomaly 

D0U3LS  PRECISION  E,TA,EA,PI 

EA  ■  DACCS((F.  +  DCOS(TA))/(  1.  OD+OO  +  E*DCOS(TA))) 
IF  (TA  .GT.  PI)  THEN 
EA  «  (2*PI)  -  EA 
ENDIF 
RETURN 
END 

SUBROUTINE  MSANMO(A,MM,MU) 

*  THIS  SUBROUTINE  CALCULATES  THE  MEAN  MOTION 

DOUBLE  PRECISION  A, MM, MU 

MM  -  DSQRT(MU/(A**3) ) 

RETURN- 

END 


SUBROUTINE  MEANAN’C EA , E ,MA) 

*  THIS  SUBROUTINE  CALCULATES  THE  MEAN  Anomaly 

DOUBLE  PRECISION  EA,E,MA 

MA  ■  EA  -  E*DSIN(EA) 

RETURN 

END 


SUBROUTINE  TFLGHT( MM , MA , TF ) 

*  THIS  SUBROUTINE  CALCULATES  THE  TIME  OF  FLIGHT 

DOUBLE  PRECISION  MM,MA,TF 

T?  =  ( 1/MM)*MA 
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*  CALCULATE  UNPERTURBED  ORBIT 

i>i>rtrtYn'tmnKn>miihnnnM'r>Vn>Vn>Vn>iyrtVnn>rift#!VtKtrtrtiVi‘f*t<*i'h'i*#)>*##rft*'**ifc#,,V*ifc*** 


SUBROUTINE  UNPRET( DT, PER , AL , LAN , AP , I , RX , R J , RK , R , 

+  vi.vj.vk.v.mu.pi^.a.e.s.ta.p.mm.ma.ea, 

t  TF ,  T ,  SUM ,  R I  RAY ,  R  JRAY ,  RKRAUf ,  RARAY , TARAY ,  AINRAY ,  APRAY , TI MRAY , 

+  IT) 

*  THIS  SUBROUTINE  CALCULATE  THE  UNPERTURBED  ORBIT 

*  THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 

*  NEWELT  «  CALCULATE  NEW  ELEMENTS  AFTER  TIME  STEP 

*  NEVPOS  *  CALCULATE  NEW  POSITION  AFTER  TIME  STEP 

*  NEWEL  «•  CALCULATE  NEW  VELOCITY  AFTER  TIME  STEP 

*  STORE  -  STORES  POSITION  IS  ARRAYS 

DOUBLE  PRECISION  T,DT,PER,AL,LAN,AP,I ,RI,RJ,RK,R,VI ,VJ,VK,V, 

+  MU, PI .H.A.E.N.TA.P.MM.MA.EA, IT, TT 

DIMENSION  RARAYC500) ,TARAY(500) ,RIRAY(500) ,RJRAY(500) , 

+  RKRAY(500),AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

*  SET  TRUE  ANOMALY  TO  NEGATIVE  SO  LOOP  CAN  BE  EXECUTED 
IF  (TA  .GT.  6.21)  THEN 

TA  «  TA  -  (2*?I) 

ESDI? 

»>  CONTINUE  AROUND  ORBIT  TILL  CLOSE  TO  PERIGEE 
230  IF  ((TA  .LE.  6.21)  .AND.  (T  .  LE.  PER))  THEN 

*  Iacr«Ritnt  TRUE  TIME 

IT  «  TT  +  DT 

CALL  NEVELTC MM ,  MA ,  E ,  E  A ,  TA ,  TF ,  DT ,  P I ,  PER ) 

CALL  N?OS(RI ,RJ,RK,R,LAN,AP, I ,  TA,A,E) 

CALL  NVEL( E , P ,TA , LAN , AP , I , VI , VJ , VK , V , MU) 

*  INCREMENT  STEP  COUNTER  AND  STORE  VALUES 

NUM  *  NUN  +  1 

CALL  STORE ( RI ,R J , RK , R ,TA , RIRAY , R JRAY , RKRAY , 

+  RARAY ,TARAY , NUM , I , AP , AINRAY , APRAY , 

+  TT,TIMRAY) 

*  INCREMENT  TIME  STEP  COUNTER 

T*  T  +  DT 

GOTO  230 
ENDIF 
RETURN 
END 


*  CALCULATE  THE  UNPERTURBED  NEW  ELEMENTS 

Vii>Wnh!lViVinn:ih:!>iWriWrfriWrt>WnVVnV;hhtiHtVi(m:itiWriViVi(innWn>iHiHt*i(inntiV*im 


ORB  1 1700 

CRB 11710 

0RB1172O 

0RBU730 

ORB 11 740 

0RII1750 

0RB11760 

0RB11770 

0RB11780 

ORB 11790 

ORB 11800 

ORB11S10 

ORB 11820 

0RB11630 

ORB 11840 

ORBUBSO 

0RB11860 

0RBU870 

0RB11880 

ORB 11890 

ORB 11900 

ORB 119 10 

0RB11920 

ORB 11930 

ORB 11940 

0RBU9S0 

0RB11960 

0RB11970 

0RB11980 

0RB11990 

ORB  12000 

ORB 120 10 

ORB 12020 

ORB12030 

0RB12040 

0RB12050 

0RB12060 

0RB12070 

ORB12080 

ORB12090 

0RB12100 

0RB12110 

ORB12120 

0RB12130 

0RB12140 

0RB12150 

0RB12160 

0RB12170 

0RB12180 

0RB12190 

ORB12200 

0RB12210 

ORB12220 

ORB1223C 

0RB12240 

ORB  3.2250 


4" 


8V2KCITISE  NEVE LT( MM , MA , £ , E A , TA , TF , DT , PI , PER ) 

*  THIS  SUBROUTINE  CALCULATES  THE  Unperturbed  NEW  ELEMENTS 

*  THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 

*  NSA  -  NEW  ECCENTRIC  ANOMALY 

*  NTA  -  NEV  TRUE  ANOMALY 

DOUBLE  PRECISION  MM,MA,E,EA,TA,TF,DT,PI ,PER 

*  Increment  TIME  OF  FLIGHT  AND  CHECK  IF  TF  GREATER  THAN  PERIOD 
TF  «  TF  +  DT 

IF  (TF  .GT.  PER)  THEN 
TF  «  TF  -  PER 
END  IF 

*  CALCULATE  MEAN  ANOMALY  AND  USE  TO  FIND  ECCENTRIC  Anomely  THEN  NEW 

*  TRUE  ANOMALY 
MA  «  MM*(TF) 

CALL  NEA(MA,E,EA) 

CALL  NTA(£A,E,TA,PI) 

RETURN- 

END 


*  CALCULATE  PERTURBED  ORBIT 


SUBROUTINE  PRETURC DT , PER , AL , LAN , AP , I , R I , R J , RK , R , 

+  VI ,VJ,VK,V,FR,FS,FW,MU,PI ,H,A,E,N,TA,P,MM,MA,EA, 

+  TF ,T , NUM , RIRAY , RJRAY , RKRAY , RARAY ,TARAY , AINRAY , APRAY ,TIMRAY , 

+  TT , TFE A , TFSU , TFMO , TFDR A , TD I , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP ) 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBED  ORBIT. 


* 


*  THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 

*  TFORCE  *  CALCULATE  THE  TOTAL  PERTURBING  FORCE  ON  THE  SATELLITE 

*  PNEVEL  «  CALCULATE  THE  Perturbed  NEW  ELEMENTS 

*  NPOS  «  NEW  POSITION  AFTER  TIME  STEP 

*  NVEL  «  NEW  VELOCITY  AFTER  TIME  STEP 

*  PERIOD  -  PERIOD  OF  PERTURBED  ORBIT 

*  STORE  «  STORE  POSITION  AND  ELEMENTS  IN  ARRAYS  FOR  PLOTTING 

DOUBLE  PRECISION  T,DT,PER,AL,LAN,AP,I ,RI ,RJ,RK,R,VI ,VJ,VK,V, 

+  FR,FS,FV,MU,PI ,H,A,E,N,TA,P,MM,MA,EA,TF,TT, 

+  DI ,DA,DE,DMM,DMA,DLAN,DH,DAP,EI ,EJ,EK,HI ,HJ,LP,M, 

+  DVR , DVS , DVW , DVI , DVJ , D VK 

DIMENSION  RARAYC500) ,TARAY(500) ,RIRAY(500) ,RJRAY(500) , 

+  RKRAY(SOO) ,AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

*  SET  MEAN  RADIUS  OF  EARTH 

RE  *  6400.0 

DT  *  PER/50 

T  =  DT 

IF  (TA  .GT.  6.21)  THEN 
TA  «  TA  -  ( 2*PI ) 
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sndif 

IF  (7?  . G£.  ?th\  THEM 
7F  =  IF  -  FEF 
ESDI? 


CONTINUE  Around  ORBIT  FOR  ONE  PERIOD 
IF  ((TF  .LT.  PER)  .AND.  (T  .  LT.  PER))  THEN 

INCREMENT  TRUE  TIME 

TT  *  TT  +  DT 

CALL  TFORCEUL.LAN.AP.I.RIUJ.RK.R.VI.VJ.VK.V, 

+  TT,FR,FS ,FW,MU,PI , 

+  FEA , FSU , FMO , FDRA , FOR , 

+  EI.EJ.EK.E.AJ.LP.TA.PER.EA.MA.TF.P, 

+  MM,N,H,HI,HJ,DT) 

CALL  PNEWELtFR,F5,FW,H,R,A,E,>\TA,DT,I,LANtAL, 

+  AP,P,MM,MA,EA,TF ,T,MU,PI , 

+  DI,DA,DE,DMM,DMA,DUN,DH,DAP) 

CALL  NPOS( RI , R J , RK , R , LAN , AP , I ,  TA,A,E) 

CALL  NVELCE.P.TA.LAN.APJ^VI.VJ.VK.V.MU) 

CALCULATE  NEW  PERIOD  AND  RESET  TIME  STEP  AND  TIME  COUNTER 
IF  NOT  AT  END  OF  ORBIT 
IF  (T  . LT.  (PER-DT))  THEN 
CALL  PERIODU.PER.PI.MU) 

DT  »  PER/50 
T  *  TF 
END  IF 

INCREMENT  STEP  COUNTER 

NUN  »  NUM  +  1 

CALL  STORE! RI , RJ , RK , R ,TA , RIRAY , RJRAY , RKRAY , 
t  RARAY , TARAY , NUM , I , AP , AINRAY , APRAY , 

+  TT.TIMRAY) 


TOTAL  ELEMENT  CHANGES 
TDI  «  TDI  +  SNGL(ABSCDI)) 

TDA  *  TDA  +  SNGL(ABSCDA)) 

TDE  *  TDE  +  SNGLCABS(DE)) 

TDMM  *  TDMM  +  SNGL(ABS(DMM)) 
TDMA  »  TDMA  +  SNGL(ABS(DMA)) 
TDLAN  «  TDLAN  +  SNGL(ABSf DUN)) 
TDH  ■  TDH  +  SNGL(ABS(DH)‘ 

TDAP  »  TDAP  +  SNGL(ABS(DAP)) 

TFEA  »  TFEA  +  FEA 

TFSU  «  TFSU  +  FSU 

TFMO  «  TFMO  +  FMO 

TFDRA  «  TFDRA  +  FDRA 


CHECK  FOR  IMPACT 
IF  (R  .LE.  RE)  THEN 

PRINT*, ' SATELLITE  WILL  IMPACT  THE  EARTH!!  ' 
T  »  PER 

END  IF 


INCREMENT  TIME  COUNTER 
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T  =  T  +  DT 
GOTO  240 
END  IF 
RETURN* 

END 


*  CALCULATE  THE  PERTURBING  FORCES 


SUBROUTINE  TFORCE(AL,LAN,AP,I,RI ,RJ,RK,R,VI ,VJ,VK,V,TT, 

+  FR , FS ,FV , MU , PI , FEA , FSU , FMO i FDRA , FOR , 

+  El ,EJ,EK,E,A,T,LP,TA,PER,EA,MA,TF,P, 

+  MM,N,H,HI ,HJ,DT) 

*  THIS  SUBROUTINE  SUHS  ALL  THE  PERTURBING  FORCES  FOR  THE  TOTAL 

*  PERTURBING  FORCE. 

*  THE  FOLLOWING  SUBROUTINES  WERE  CALLED: 

*  OBERT  «  OBLATENESS  OF  THE  EARTH 

*  FSUN  »  GRAVITATIONAL  Ac trace ion  OF  THE  SUN 

*  FMOON  »  GRAVITATIONAL  Attraction  OF  THE  HOON 

*  FDRAG  ■  DRAG  FORCES 

DOUBLE  PRECISION  FER.FES.FEVJSR.FSSJSW.FMR.FMS.FMV.MU.PI, 

+  FDR,FDS ,FDW,FR,FS,FW,RI ,RJ,RK,R,AL,I ,TT,LAN,AP>VI ,VJ,VK,V, 
+  EI.EJ.EK.E.AJ.LP.TA.PER.EA.MA.TF.P, 

+  MM,N,H,HI,HJ,DT 


CALL  OBEARTC RI , R J , RK , R , AL , I , FER, FES ,FEW , HU) 

CALL  FSUNCTT , RI , RJ , RK , R , FSR , FSS , FSW , PI ) 

CALL  FMOON'C TT , RI , RJ , RK , R , FNR  rMS,FMV,PI) 

CALL  FDRAG( RI , R J , RK , R , VI , VJ  •  ■  /,LAN,AP,I,FDR,FDS,FDW, 

+  El  ,EJ,EK,E  *,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 

+  MM,N,H,HI,HJ,DT) 


*  SUM  VECTOR  FORCES 

FR  *  FER  +  FSR  +  FMR  +  FDR 
FS  *  FES  +  FSS  +  FMS  +  FDS 
FW  *  FEW  +  FSW  +  FMW  +  FDW 

*  CALCULATE  TOTAL  FORCE  FROM  EACH,  AND  TOTAL  OF  ALL 
FEA  «  SNGL(SQRT((F£R**2)+(FES**2)+(FEW**2))) 

FSU  «  SNGL(SQRT( (FSR**2)+(FSS**2)+(FSW**2) ) ) 

FMO  «  SNGL(SQRT((FMR**2)+(FMS**2)+(FMW**2) ) ) 

FDRA  »  SNGL(SQRT((FDR**2)+(FD5**2)'KFDW**2))) 

FOR  *  SNGL(SQRT( (FR**2)+(FS**2)+(FW**2) ) ) 


RETURN 

END 

irkhiiltltirkirklrtrlrk'k'kirktr* 


SUBROUTINE  OBEART( RI , R J , RK , R , AL , I , FER , FES , FEW , MU) 

*  THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  THE 

*  OBLIQUENESS  OF  THE  EARTH. 
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SUBROUTINE  FSUN(TT,RI,RJ,RK,R,FSR,FSS, FSW , PI ) 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  THE  SUN 

THE  FOLLGVING  SUBROUTINES  ARE  CALLED: 

SUNPOS  *  SUNS  POSITION  ORBITING  AROUND  EARTH 
HEVBOD  *  PERTURBING  FORCE  FROM  A  Httvtnly  BODY 

DOUBLE  PRECISION  FSR,FSS,FSW,PI, 

+  RSI.RSJ.RSK.SLAN.SI.SAL.SMU.TT.RI.RJ.RK.R.RS 

SUNS  GRAVITATIONAL  PARAMETER 
SMU  «  1. 3271544D+11 

CALL  SUNPOS ( TT , RS I , RS J , RSK , RS , SLAN , S I , S AL , P I ) 


RETURN- 

END 


DOUBLE  PRECISION  JS.RE.FER.FES.FEW.RI.RJ.RK.R.AL.I.MU.M  CRB13940 

ORBI39SO 

J2  *  1.P82364D-W  0RB13960 

RE  *  6.37S2E+03  ORB13970 

0RB13980 

PER  *  ((-3.QD4-00*HU*J2*(R£**2))/(2.0IH00ft(Rft*4)))*  0RB13990 

+  (1.0D+00  -  (3.  0D4-00*  ((DSIN(  I))**2)  *  ((DSIN(AL))**2)))  0RBI4000 

FES  ■  (  ( « 3.  GD+00*MU*J2*( RE**2) )/(R**4) )*  QRB14010 

+  (C(DSIX(I))**2)*  (DSIN(AL))*(DCOS(AL)))  0RB14020 

FEW  *  ((-3.0D+00*MU*J2*(RE**2))/(R**4))*  ORB14030 

+  ( DS IN( I )*DCOS( I )*DSIN( AL) )  ORB  14040 

RETURN  ORB 14050 

ENO  ORB 14060 

.  0RB14070 

0RB14090 

SUBROUTINE  FSUN(TT,RIfRJ,RK,R,FSR,FSS,FSV,PI)  0RB14100 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  THE  SUN  ORB14110 

0RB14120 

THE  FOLLGVING  SUBROUTINES  ARE  CALLED:  ORB 14 130 

SUNPOS  *  SUNS  POSITION  ORBITING  AROUND  EARTH  ORB 14 140 

HEVBOD  *  PERTURBING  FORCE  FROM  A  Httvtnly  BODY  0RB14150 

ORB 14 160 

DOUBLE  PRECISION  FSR,FSS,FSW,PI,  0RB14170 

K  RSI.RSJf.RSK.SLAN.SI.SAL.SMU.TT.RI.RJ.RK.R.RS  ORB14180 

ORB14190 

SUNS  GRAVITATIONAL  PARAMETER  ORB 14200 

SMU  «  1.3271544D+U  0RB14210 

CALL  SUNPOSCTT.RSI ,RSJ,RSK,RS,SLAN,SI ,SAL,PI)  0RB14220 

CALL  KEVBOD(RI ,RJ,RK,R,RSI,RSJIRSK,RS,SLAN,SAL,SI,SMU,FSR,FSS,FSW)ORB14230 
RETURN  ORB 14240 

END  ORB 14250 

.  ORB  14260 

ORB14280 

SUBROUTINE  FMOON(TT,RI , R J , RK , R , FMR , FMS ,FMV , PI )  ORB 14290 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  Tht  MOON  0RB14300 

0RB14310 

THE  FOLLOWING  SUBROUTINE  ARE  CALLED:  ORB 14320 

MONPOS  *  MOONS  POSITION  ORBITING  AROUND  THE  EARTH  ORB14330 

HEVBOD  *  PERTURBING  FORCE  FROM  A  HEAVENLY  BODY  0RB14340 


*  THE  FOLLOWING  SUBROUTINE  ARE  CALLED:  ORB 14320 

*  MONPOS  *  MOONS  POSITION  ORBITING  AROUND  THE  EARTH  ORB14330 

*  HEVBOD  *  PERTURBING  FORCE  FROM  A  HEAVENLY  BODY  ORB14340 

ORB14350 

DOUBLE  PRECISION  FMR.FMS.FMV.RMI.RMJ.RMK.MLAN.MI.MAL.MMU,  ORB14360 

+  TT,RI ,RJ,RK,R,RM,PI  ORB14370 

ORB 14380 

*  MOONS  GRAVITATIONAL  PARAMETER  0RB14390 

MMU  a  4. 90287D+03  ORB 14400 

0RB14410 

CALL  MONPOS(TT,RMI ,RMJ,RMK,RM,MLAN,MI ,MAL,PI)  0RB14420 

CALL  HEVBODC RI , RJ , RK , R , RMI , RM J , RMK , RM , MLAN ,MAL,MI, MMU , FMR , FMS , FMV) ORB 14430 
RETURN  ORB 14440 

END  ORB 14450 

ORB14460 
ORB  14470 
ORB 14480 

SUBROUTINE  HEVBOD( RI , R J , RK , R , RPI , RP J , RPK , RP , LAN , AL , INC , MUP ,  ORB 14490 


+  FHR,FHS,FKW)  GRB 14500 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  A  QRB14510 

HEAVENLY  BODY.  ORB 145 20 

ORB  14530 

THE  FOLLOWING  SUBROUTINE  WAS  CALLED:  0RB14540 

IJKSSV  *  1 IJK1  SYSTEM  TO  THE  * RSW*  SYSTEM  ORB14550 

ORB 145 60 

DOUBLE  PRECISION  DOT.FHI.FHJ.FHK.RI.RJ.RK.R.RPI.RPJ.RPMP,  0RB14570 

+  LAN.AL.INC.MUPJ.J.K.IP.JP.KP.M.FKR.FHS.FKV  ORB14580 

0RB14590 

CALCULATE  UNIT  VECTOR  FOR  SATELLITE  AND  PERTURBING  BODIES  POSITIONORB 14600 


I  *  RI/R 
J  *  RJ/R 
K  *  RK/R 
IP  -  RPI/RP 
J?  *  RPJ/RP 
KP  *  RPK/RP 

CALCULATE  DOT  PRODUCT  OF  UNIT  VECTORS 
DOT  *  ((  I*IP  )+(  J*JP  )+(  K*KP  )) 

CALCULATE  FORCES  IN  THE  * IJK'  SYSTEM 
FHI  *  (MUP/(RP**2) )*(K/RP)*( 3,  0D+00*D0T*( IP) -( I)) 
FHJ  *  (MUP/(RP*>>2))*(R/RP)*(3.0D+00*DOT*(JP)-(J)) 
FHS  *  (MUP/(RP**2))*(R/RP)*(3.0D+0Q*DOT*(KP)-(K)) 

Transform  FORCES  TO  THE  RSW  SYSTEM 

CALL  I JKRSW( LAN , AL , INC , FHI , FHJ , FHK , FHR , FHS , FHW ) 

RETURN 

END 


ItMrMtlrMtMrMtiti 


SUBROUTINE  SUNPOS(TT . RS I , RS J , RSK , RS , SLAN , S I , SAL , PI ) 
THIS  SUBROUTINE  CALCULATES  THE  SUNS  POSITION 

VARIABLES  USED  TO  DESCRIBE  THE  SUNS  ORBIT: 

SI  *  SUNS  INCLINATION 

SLAN*  SUNS  Longituda  OF  ASCENDING  NODE 

SAP  *  SUNS  ARGUMENT  OF  PERIGEE 

RS  =  SUNS  ORBITAL  RADIUS 

STA  *  SUNS  TRUE  ANOMALY 

SAL  *  SUNS  ARGUMENT  OF  LONGITUDE 

DOUBLE  PRECISION  SLAN.SI ,SAL,RS , STA, SAP, TT, RSI ,RSK, 
+  RSJ,RSP,RSQ,RSV,PI 

SI  =  4. 09279709D-01 
SLAN  *  O.OD+OO 
SAP  *  O.OD+OO 
RS  =  1. 4959965D+08 

STA  *  ((2. 0*PI)/(365. 0  *  86400.0)  *  TT) 

SAL  =  STA  +  SAP 

CALCULATE  SUNS  POSITION  IN  ' PQW '  SYSTEM 
RS?  =  RS*DCOS(STA)  • 
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RSQ  «  RS*DSIN(STA) 

RSV  «  O.ODtCQ 

*  TRANSFORM  FOSITION  TO  'IJK'  SYSTEM 

CALL  PQVI JK( SLAN , S AP , S I , RSP , RSQ , RSW ,  RS  I , RS J , RSX ) 

RETURN 

END 


*********************i‘t********'******  **********  ft***  ft***  ********** 

SUBROUTINE  MONPOS(TT, RMI , RMJ , RMK , RM, MLAN, MI ,MAL# PI ) 

*  THIS  SUBROUTINE  CALCUUTES  THE  MOONS  POSITION 

*  VARIABLES  USED  TO  DESCRIBE  THE  SUNS  ORBIT: 

*  MI  »  MOONS  INCLINATION 

*  MLAN-  MOONS  Longitude  OF  ASCENDING  NODE 

*  MAP  *  MOONS  ARGUMENT  OF  PERIGEE 

*  RM  »  MOONS  ORBITAL  RADIUS 

*  MTA  *  MOONS  TRUE  ANOMALY 

*  MAL  *  MOONS  ARGUMENT  OF  LONGITUDE 

DOUBLE  PRECISION  8  MLAN, MAL, RM,TM, MTA, RMP, RMQ, RMV, 

+  RMI, RMJ, RMK, L  .-,PI 

MI  «  4.99164166D-01 
RM  -  3. 844D+05 
MLAN  *  0.  0 

MTA  «  ((2. 0*PI)/(27. 3  *  3600)  *  IT) 

MAP  *  O.OD+OO 
MAL  *  MTA 

*  CALCULATE  MOON  POSITION  IN  'PQV'  SYSTEM 
RMP  «  RM*DCOS(MTA) 

RMQ  «  RM*DSIN(MTA) 

RMV  «  0 

*  TRANSFORM  POSITION  TO  'IJK'  SYSTEM 

CALL  PQVI JK( MLAN , MAP , MI , RMP , RMQ , RMV , RMI , RMJ , RMK) 

RETURN 

END 

***"t'  'ft*********************************************************** 


SUBROUTINE  FDRAG(RI ,RJ,RK,R,VI ,VJ,VK,V,LAN,AP,I ,FDR,FDS,FDW, 

+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 

+  MM,K,H,HI,HJ,DT) 

*  THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  DRAG 


Vf 

* 

Vr 

* 

* 


THE  FOLLOWING  VARIABLES  ARE  USED  TO  MODEL  THE  ATMOSPHERE: 

RE  *  RADIUS  OF  EARTH 

M  =  MASS  OF  SATELLITE 

AR  *  FRONTAL  SURFACE  AREA  OF  SATELLITE 

Z  =  ALTITUDE  OF  SATELLITE 

K  =  EXPONENTIAL  DECAY  FACTOR 

DENO  =  NORMAL  DENSITY 

CD  COEFFICIENT  OF  DRAG 


0RB15060 

0RB15070 

0RBIS080 

ORB 15090 

0RB15100 

ORB IS 110 

ORB15120 

0RB15130 

0RB15140 

0RB15150 

0RB15160 

0RBI5170 

ORB 13 180 

0RB1S190 

ORB 15200 

0RB15210 

ORB 15 220 

0RB15230 

ORB 13240 

0RB1S250 

ORB 15260 

0RB15270 

0RJ13280 

0RB15290 

ORB 15300 

0RB15310 

0RB1S320 

0RB1S330 

ORB 13340 

ORB15350 

0RB15360 

0RB15370 

ORB15380 

0RB15390 

ORB 15400 

ORB15410 

0RB15420 

ORB 15430 

ORB 15440 

ORB 15450 

ORB 15460 

ORB 15470 

0RB15480 

0RB15490 

ORB 15500 

0RB15510 

0RB15520 

0RB15530 

ORB 15540 

0RB15550 

0RB15560 

0RB15570 

0RB15580 

0RB15590 

ORB15600 

ORB 15610 
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DOUBLE  PRECISION 

* 

+ 


MAG,M,K,FDR,FDS,FDV,RE,AR,Z,DENO,CD,DEN. 

FDJ.FDKJDI.RI.RJ.RK.VI.VJ.VK.V.LAN.AP.I.R, 

El  ,EJ,EK,E,A,T,LP,TA,FF,R,EA,MA,AL,TF,P,PI  ,MU, 
MM ,N,K,HI ,HJ,DT, DVR , DVS ,DW , DVI , DVJ , DYK 


RE  «  6.378145D+03 
M  »  1.  0D+02 
AR  *  2,00+01 
Z  ■  R  -  RE 

*  DEPENDING  ON  ALTITUDE  SET  ATMOSPHERE  VARIABLES 
I?  (Z.LE.  1.5D+02)  THEN 

K  -  4. 74D-02 
DENO  «  1. 225D+00 
CD  «  1.  OD+OO 

ELSEIF  (Z.LE.5.5D+02)  THEN 
K  «  3.4614D-02 
DENO  -  1.79846D-01 
CD  «  2.  OD+OO 
ELSE 

K  *  2.21698D-3 
DENO  «  1.015484D-07 
CD  «  2. OD+OO 
END  IF 

*  CALCULATE  ATMOSPHERIC  DENSITY 
DEN  *  DENO  *  DEXP(-K*Z) 

*  CALCULATE  MAGNITUDE  OF  DRAG  FORCE  AND  LIMIT  IT  TO  1. OE-20 
MAG  «  -C0.5D+00)*CD*AR*DEN*V*(l.QD-03)/M 

I?  (ABS(MAG)  .LT.  1.0D-20)  THEN 
MAG  «  -1.0D-20 
ENDIF 

*  GIVE  DRAG  FORCE  A  Direction  OF  MINUS  THE  VELOCITY 
FDR  »  0.  0 

FDS  «  MAG  *  V 
FDW  «  0.  0 
RETURN 
END 


*  CALCULATE  PERTURBED  NEW  ELEMENTS 


SUBROUTINE  PNEWEL(FR,FS,FV,H,R,A,E,N,TA,DT,I,LAN,AL,AP,P, 

+  MM , MA , EA ,TF ,T , MU , PI , DI , DA , DE , DMM , DMA , DLAN , DH , DAP) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  ELEMENTS  FROM  THE  PREVIOUS 
ELEMENTS  ADDED  TO  THE  RATES  OF  CHANGE  FOR  ONE  STEP 

THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

RATE  =  CALCULATES  RATES  OF  CHANGE  OF  ORBITAL  ELEMENTS 
NANGMO  =  NEW  ANGULAR  MOMENTUM  (NEVH) 

NSMA  =  NEW  SEMI -MAJOR  AXIS  (NEWA) 

NECC  =  NEW  ECCENTRICITY  (NEWE) 


0RB15620 
0RB15630 
0RB13640 
0RB15650 
ORB 15660 
0RB15670 
0RB15680 
0RB15690 
0RB15700 
ORB 15 7 10 
ORB15720 
ORB15730 
0RB15740 
ORB15750 
0RB15760 
0RB15770 
ORB15780 
0RB15790 
ORB1S800 
0RB15810 
0RB15820 
0RB15S30 
ORB 15840 
ORB15850 
0RB15860 
ORB15870 
0RB15860 
0RB15890 
ORB 15900 
0RB15910 
0RB15920 
ORB 15930 
ORB 15940 
ORB15950 
0RB15960 
0RB15970 
0RB15980 
0RB15990 
ORB 16000 
ORB 160 10 
0RB16020 
ORB 16030 
ORB16040 
ORB 16050 
ORB 16060 
ORB16070 
ORB 16080 
ORB16090 
0RB16100 
0RB16110 
0RB16120 
0RB16130 
0RB16140 
0RB16150 
0RB16160 
0RB16170 
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*  S'IS'CL  «  NEW  INCLINATION  (NEtfl) 

*  NASNOD  *  NEW  LONGITUDE  OF  ASCENDING  NODE  (NEVLAN) 

*  NASPSR  *  NEW  ARGUMENT  OF  PERIGEE  (  NEVA?) 

*  NNNNO  »  NEW  MEAN  MOTION  (NEWMM) 

*  MSANMO  «  MEAN  MOTION  (MM) 

*  NMMAN  «  NEW  MEAN  ANOMALY  (NEVMA) 

*  NEA  «  NEW  ECCENTRIC  ANOMALY  (EA) 

*  NTA  «  NEW  TRUE  ANOMALY  (TA) 

*  TFLGKT  »  TIME  OF  FLIGHT  (TF) 

DOUBLE  PRECISION  FR,FS,FV,DMM,H,R,A,E,N,TA,DT,I,LAN,AL,AP,P, 
+  MM ,MA , EA ,TF ,T, MU , PI , DA , DH , DE , DI , D LAN , DAP , DMA , 

+  NEWH , NEVA , NEVE ,  NEVI , NEVLAN , NEWAP , NEVMM 

*  INCREMENT  TIME  OF  FLIGHT  BY  ONE  TIME  STEP  AND  CALCULATE  RATES 
TF  *  TF  +  DT 

CALL  RATESCDH.DA.DE.DI.DLAN.DAP.DMM^MA.E.MMjR.A.FRJS.FV, 

+  TA,AL,H,P,T,MU, I) 

*  CALCULATE  NEW  ELEMENTS 
CALL  NANGMO( H , DT , DH , NEWH ) 

CALL  NSMAC  A , DT , DA , NEVA ) 

CALL  NECC( E , DT , DE , NEVE ) 

CALL  NINCLC I , DT , DI , NEW I ) 

CALL  NASNODC  LAN , DT , DLAN , NEVLAN) 

CALL  NARPER(AP,DT, DAP, NEWAP) 

*  SET  ELEMENTS  TO  NEW  ELEMENTS 
A  *  NEVA 

E  *  NEVE 
I  *  NEW I 
LAN  *  NEVLAN 
AP  *  NEWAP 
P  «  A  *  (1  -  E**2) 

*  MOVE  THE  SATELLITE  ONE  TIME  STEP 
CALL  MEANMO(A,MM,MU) 

CALL  NMNAN( MA , MM , DT ,TF , DMA , PI ) 

CALL  NEA(MA,E,EA) 

CALL  NTA(EA,E,TA,PI) 

CALL  TFLGHT( MM , MA ,TF) 

AL  »  TA  +  AP 

RETURN 

END 


*  CALCULATE  THE  RATES  OF  CHANGE  OF  THE  ORBITAL  Elements 


SUBROUTINE  RATESCDH.DA.DE.DI.DLAN.DAP.DMM^MA^MM.R.A.FR^S.FV, 
+  TA,ALkH,P,T,MU,I) 

THIS  SUBROUTINE  Calls  THE  FOLLOWING  SUBROUTINES  TO  CALCULATE  THE 
TIME  RATE-OF-  CHANGE  OF  THE  ORBITAL  ELEMENTS: 

RSMAX  =  RATE-OF-CHANGE  OF  THE  SEMI-MAJOR  AXIS  (DA) 

RECC  =  RATE-OF-CHANGE  OF  THE  ECCENTRICITY  (DE) 

RINC  =  RATS-Gr -CHANGE  OF  THE  INCLINATION  (DI) 


0RB16180 
0RB16190 
QRB16200 
0RBI6210 
0RB16220 
0RB16230 
0RB16240 
ORB 16250 
0RB16260 
0RB16270 
ORB 16280 
0RB16290 
ORB 16300 
ORB 163 10 
0RB16320 
0RB16330 
ORB 16340 
0RB16350 
ORB16360 
0RB16370 
0RB16380 
0RB16390 
ORB 16400 
0RB16410 
ORB 16420 
ORB 16430 
ORB 16440 
0RB16450 
ORB 16460 
0RB16470 
ORB 16480 
ORB 16490 
ORB 16500 
0RB16510 
ORB 16520 
0RB16530 
ORB 16540 
ORB 16550 
ORB 16560 
0RB16570 
ORB165BO 
0RB16590 
ORB 16600 
0RB16610 
ORB16620 
ORB 16630 
ORB 16640 
0RB16650 
ORB16660 
0RB16670 
ORB16660 
0RB16690 
0RB16700 
0RB16710 
ORB16720 
ORB16730 


*  RLAN  »  RATE-QF-CHANGE  OF  THE  Lonjitud*  OF  THE  ASCENDING  NODE 

*  (DLAN) 

*  SAP  *  RATE-OF-CHANGE  OF  THE  ARGUMENT  OF  PERIGEE  (DAP) 

*  EMM  *  RATE-QF- Chang*  OF  THE  MEAN  MOTION  (DMM) 

*  RMA  »  RATE-OF-CHANGE  OF  THE  MEAN  ANOMALY  (DMA) 

*  RANGMO  »  RATE-OF-CHANGE  OF  THE  ANGULAR  MOMENTUM  (DH) 

DOUBLE  PRECISION  DH,DA,DE,DI ,DLAN,DAP,DMM,DMA,E,MM,R,A,FR,FS,FW, 
+  TA,AL,H,P,T,MU,I 

CALL  R5MAX( E,MM , R , A , FR , FS , DA ,TA) 

CALL  RF.CC(E,MM,R,A,FR,FS,TA,DE) 

CALL  RINC(E,MM,R,A,FV,AL,DI) 

CALL  RLAN(E,MM,R,A,I ,FV,AL,DLAN) 

CALL  RAP(E,MM,R,A,I,H,P,AL,TA,FR,FSIFV,DAP) 

CALL  RMM( MM , A , DMM , DA , MU) 

CALL  RMA(E,MM,R,A,TA,DMM,FR,FSIDMA,T) 

CALL  RANGMO(R , FS , FV , DH) 

RETURN 

END 

ft**-*frfrft-ftftVn>*Tfc*ft******tff#***fr*********************6************fr* 

SUBROUTINE  RANGMO ( R , FS , FV , DH) 

*  THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE 

*  ANGULAR  MOMENTUM 

DOUBLE  PRECISION  FS,FV,DHtf,DKS,DK,R 


DHV  «  R  *  FS 
DHS  »  R  *  FV 

DH  «  DSQRT((DHV»'Wr2)  +  (DHS**2)) 

RETURN 

END 


SUBROUTINE  RSMAX( E , MM , R , A , FR , FS , DA ,TA) 

■>  THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  SEMI -MAJOR 

*  AXIS 

DOUBLE  PRECISION  DA,FR,FS ,E,MM,R>A,TA,ET 

*  TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (E.GT.0.9)  THEN 

ET  *  0.  9 
ELSE 

ET  *  E 
ENDIF 

DA  =  ((2. 0D+00*E  *DSIN(TA))/(MM  *DSQRT(1.  0D+00-(ET**2))))*FR  + 
+  ((2.  0D+00*A*DSQRT(1.  0D+00-(E  **2) ) )/(MM  *R))*FS 

RETURN 
END 


0RB16740 
0RB16750 
ORB16760 
0RBI6770 
0RB16780 
0RB16790 
ORB168O0 
0RB16810 
ORB16820 
ORB 16830 
ORB 16840 
ORB168SO 
ORB 16860 
0RB16870 
ORB 16880 
ORB16890 
ORB16900 
0RB16910 
0RB16920 
0RB16930 
ORB 16940 
ORB 16950 
ORB 16960 
0RB16970 
ORB 16980 
0RB16990 
0RB17000 
ORB17010 
ORB17020 
0RB17030 
ORB 17040 
ORB17050 
0RB17060 
0RB17070 
0RB17080 
ORB 17090 
ORB 17 100 
0RB17110 
0RB17120 
ORB  17 130 
0RB17140 
0RB17150 
0RB17160 
0RB17170 
ORB 17 180 
ORB17190 
ORB17200 
0RB17210 
ORB17220 
0RB17230 
0RB17240 
0RB17250 
0RB17260 
0RB17270 
0RB17280 
0RB17290 
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* 

* 


* 


SUBROUTINE  RECC( E,MM f R , A , FR , FS ,TA , DE) 

*  THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  ECCENTRICITY 

DOUBLE  PRECISION  DE , FR , FS , E , MM , R , A ,TA , ET 

*  TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  2ERO 
IF  (E.LT.0. 1)  THEN 

ET  *  0.  1 
ELSE 
ET  *  E 
ENDIF 

DE  -  «DSQRT(  1.  OD+OO  -  (E  **2))*SIN(TA))/(MM  *A))*FR  + 

+  (CDSQRTC 1.  OD+OO  -  (E  **2)))/(MM  *ET*(A**2)))* 

+  ((A**2)*(l. OD+OO  •  (E  **2))/(R)  -  (R))*FS 

RETURN 
END 


SUBROUTINE  RLAN( E , MM , R , A , I , FV , AL , DLAN) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  LONGITUDE 
OF  THE  ASCENDING  NODE 

DOUBLE  PRECISION  DLAN,FV,E,MM,R,A,I ,AL,ET,IT 

TRAP  (E)  AND  (I)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (E.GT.0.9)  THEN 
ET  «  0.  9 
ELSE 
ET  «  E 
ENDIF 

IF  ( I.  LT.  0.01745)  THEN 
IT  *  0.01745 
ELSE 

IT  *  I 
ENDIF 

DLAN  «  (R"FW*DSIN(AL))/(MM  *(A**2)*DSQRT( 1.  OD+OO  -  (ET**2))* 

+  DSIN(IT)) 

RETURN 

END 

SUBROUTINE  RAPCE.MM.R.A.I.H.P.AL.TA.FR.FS.FV.DAP) 

*  THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  ARGUMENT 

*  OF  PERIGEE 

DOUBLE  PRECISION  DAPR,DAPS,DAPV,DAP,FR,FS,FV,E,MM,R,I,H,P,AL,TA, 
+  ET.A.IT 

*  TRAP  (E)  AND  (I)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (I.  LT.  0.01745)  THEN 

IT  =  0.01745 
ELSE 

IT  =  I 
ENDIF 


0RB17300 

ORB 173 10 

ORB  17320 

0RB17330 

0RB17340 

0RB17350 

0RB17360 

0RB17370 

ORB  17  380 

ORB17390 

0RB17400 

0RB17410 

ORB17420 

ORB17430 

ORB17440 

0RB17450 

0RB17460 

0RB17470 

ORB 17480 

ORB 17490 

ORB17500 

ORB17510 

0RB17520 

0RB17530 

0RB17540 

0RB17550 

ORB 17560 

0RB17570 

OR317580 

ORB17590 

ORB17600 

0RB17610 

0RB17620 

ORB17630 

0RB17640 

ORB17650 

ORB17660 

0RB17670 

ORB17680 

ORB 17690 

ORB17700 

0RB17710 

0RB17720 

0RB17730 

ORB17740 

0RB17750 

ORB17760 

0RB17770 

0RB17780 

ORB17790 

ORB17800 

ORB17810 

0RB17820 

0RB17830 

ORB 17 840 

0RD17850 
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IF  (E.GT.0.3)  THEN 
ET  =  0.  9 

ELSEIF  CE.LT.0. 1)  THEN 
ET  »  0.  1 
ELSE 
FT  *  E 
ENDIF 

DAPR  -  (-DSQRT( 1.0+00  -  (E  **2))*DCOS(TA))/(MM  *A*ET)  *  FR 
DAPS  -  (P/(£T*H))*(DSIN(TA))* 

+  (1. OD+OO  +  1.  0D+00/(l.  OD+OO  +  ET'DCOS(TA)))  *FS 

DAPV  -  (-R*(  1.  0D+00/DTAN(  IT)),VDSIN(AL))/ 

+  (MM  *(A'V*2)*DSQRT(  1.  OD+OO  -  (ET>>*2)))*FV 
DAP  *  DAPR  +  DAPS  +  DAPV 
RETURN 
END 


SUBROUTINE  RINC( E , MM , R , A , FV , AL , DI ) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  INCLINATION 

DOUBLE  PRECISION  DI ,FV,E,MM,R,A,AL,ET 

TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  2ER0 
IF  (E.GT.0.9)  THEN 
ET  «  0.  9 
ELSE 

ET  »  E 
ENDIF 

DI  =  ( R*FW*DCOS( AL) ) / ( MM  *(A**2)*DSQRT( 1. OD+OO  -  (ET**2))) 

RETURN 

L..D 


SUBROUTINE  RMM( MM , A , DMM , DA , MU ) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  MEAN  MOTION 

DOUBLE  PRECISION  DMM, DA, MM, A, MU 

DMM  «r(-3.0D+00*MU)/(2.0D+00*MM  *(A**4)))*  DA 

RETURN 

END 


SUBROUTINE  RMA( E , MM , R , A ,TA , DMM , FR , FS , DMA ,T) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  MEAN  Anomaly 

DOUBLE  PRECISION  DMAA,DMAB,DMAC,DMAD,DMM,FR,FS,DMA,E,MM,R,A,TA, 

+  ET,T 

TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (E.GT.0.9)  THEN 
ET  =  0.  9 

ELSEIF  (E.LT.  0.  1)  THEN 


0RB17860 

0RB17870 

ORB17880 

0RB17890 

0RB17900 

0RB1791Q 

ORB 17920 

0RB17930 

0RB17940 

0RB17950 

0RB17960’ 

0RBI7970 

ORB17980 

0RB17990 

ORB18000 

0RB18010 

0RB18020 

ORB18030 

ORB 18040 

0RB18050 

0RB18060 

0RB18070 

ORB 18080 

ORB18090 

0RB18100 

0RB18110 

0RB18120 

ORB18130 

ORB 18 140 

0RB18150 

ORB 18 160 

0RB18170 

0RB18180 

ORB 18 190 

0RB18200 

0RB18210 

0RB18220 

0RB18230 

ORB 18240 

0RB18250 

0RB18260 

0RB18270 

0RB18280 

0RB18290 

ORB18300 

0RB18310 

0RB18320 

0RB18330 

0RB18340 

0RB18350 

0RB18360 

ORB18370 

0RB18380 

ORB18390 

ORB18400 

0RB18410 


5S 


F.7  3  0.  1 
SIFT. 

FT  *  S 
END  IF 

DJIA  *  (-1.0D+00/O1M  *A))* 

+  (((2.  CD-rOOVrR)/A)  -  ((1  -  (E  **2))/ET)*DC0S(TA))  *  FR  - 

+  ( 1-(E  **2))/(MM  *A*ET)*(  1+  R/(A*(  1-(E**2))))*(SIN(TA)*FS) 

+  (T  *  DMM) 

RETURN 

END 


r  TT  TiTtTTi  TiT.T.  T.TiTYXVXj  7TT> . 


*  CALCULATE  THE  NEW  ORBITAL  ELEMENTS 


,  T,T\T\TvT»TN7NTt7jT>T,T  T;7Y7  TS* 


,  ,TiTiT»T»T,T  TIT.* 


,  A  »  *  4  A  A  J  i  A  ,  ft  i  i 


SUBROUTINE  NSMA( A , DT , DA , NEVA) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  SEMI -MAJOR  AXIS 

DOUBLE  PRECISION  DA, DT, A, NEKA 

NEVA  *  A  +  DA*DT 

RETURN 

END 


.TiT.T^TrT.T  TTTTTTTXT  T  T  T  T  TiT  TTiT  i 


SUBROUTINE  NECC( E , DT , DE , NEVE) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  ECCENTRICITY 

DOUBLE  PRECISION  DE,DT,E,NEVE 

NEVE  *  E  +  DE*DT 

RETURN- 

END 

SUBROUTINE  NINCLC I ,DT,DI ,NEVI) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  INCLINATION 

DOUBLE  PRECISION  DI,DT,I,NEWI 

NEVI  «  I  +  DI*DT 

RETURN 

END 


SUBROUTINE  NASNODC LAN , DT, DLAN , NEVLAN) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  LONGITUDE  OF  THE  ASCENDING 


D0U3LE  PRECISION  DLAN ,DT, LAN, NEW LAN 

NEVLAN  =  LAN  +  DLAN*DT 

RETURN- 

END 


0RBI8420 
0RB1B430 
ORB  18440 
0RB18450 
0RB18460 
0RB18470 
-  ORB 18480 
ORB 18490 
ORB18SOO 
0RB18510 
0RB18520 
ORB 18530 
0RB18S40 
0RB18SS0 
0RB18S60 
0RB18S70 
ORB18S80 
ORB18S90 
ORB 18600 
ORB 186 10 
ORB 18620 
0RB1S630 
0RB18640 
0RB18650 
0RB18660 
0RB18670 
0RB18680 
0RB18690 
0RB18700 
ORB18710 
0RB18720 
0RB18730 
ORB18740 
ORB18750 
ORB18760 
ORB18770 
0RB18780 
ORB18790 
ORB18800 
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SUBROUTINE  NARPER(AP , DT, DAP .NEVAP) 

THIS  SUBROUTINE  CALCULATES  THE  NEW*  ARGUMENT  OF  PERIGEE 

DOUBLE  PRECISION  DAP.DT.AP, NEVAP 

NEVAP  »  AP  +  DAWY 

RETURN 

END 


ft********************************************************'****** 

SUBROUTINE  NMNAN(MA,MM,DT,TF,DMA,PI) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  MEAN  Anomaly 

DOUBLE  PRECISION  DMM,FR,FS,DMA,DT,MA,E,R,A,TA,MM,TF,T,PI 

MA  «  MH*(TF)  +  DMA*DT 
IF  (MA  .CT.  (2*PI))  THEN 
MA  -  MA  -  (2*PI) 

ENDIF 

RETURN 

END 

**************************************************************** 

SUBROUTINE  NMNMO( MM , DMM ,DT, NEVMM) 

*  THIS  SUBROUTINE  CALCULATE  THE  NEW  MEAN  MOTION 

DOUBLE  PRECISION  DMM, DT, MM, NEVMM 


NEVMM  -  MM  +  DMM*DT 

RETURN 

END 

ftftftftftftVfftftftVrftftftftftftftftftftft ft* ********  ********************************* 

SUBROUTINE  NEA(MA,E,EA) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  ECCENTRIC  ANOMOLY  BY  USING 

*  NEWTONS  METHOD  OF  ROOT  FINDING 

DOUBLE  PRECISION  EAN,MAN,MA,E,EA,DIFF 

*  LET  (EA)  EQUAL  (MA)  FOR  INITIAL  GUESS  AT  ROOT 
EA  -  MA 

SAN  *  EA  +  (MA  -  EA  4E*DSIN(EA))/(X. 0D+00  -  E*DCOS(EA)) 

MAN  *  EAN  -  E*SXN(EAK) 

*  CHECK  DIFFERENCE  (DIFF) 

DIFF  =  ABS(MA  -MAN) 

EA  *  EAN 

*  CONTINUE  TO  IXTERATE  UNTIL  DIFFERENCE  IS  NEGLIGIBLE 
200  IF(DIFF.GT.  0.0000000001)  THEN 
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SAN  *  EA  +  (HA  -  EA  +  E*DSIN(EA))/(1. OD+OO  -  E*DCQS(SA)) 
HAS  *  HAS  -  E*DSXK(EAX) 

EA  *  EAN 

DIF?  «  A3S(MA  -  HAN) 

GOTO  200 
ESDI? 

EA  -  EAN 

RETURN 

END 


SUBROUTINE  NTA(EA,E,TA,PI) 

*  THIS  SUBROUTINE  CALCUUTES  THE  NEW  TRUE  Anomaly 
DOUBLE  PRECISION  EA,E,TA,PI 

TA  «  DACOS((E  -  BCOS(EA))/(E*DCOS(EA)  -  1.  OD+OO)) 

IF  (EA. GT. PI)  THEN 
TA  ■  (2*P1)  -  TA 
END  IF 
RETURN 
END 

SUBROUTINE  NANGMO( H . DT , DH , NEWH ) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  ANGULAR  MOMENTUM 

DOUBLE  PRECISION  DH,DT,K,NEVK 

NEWH  »  K  +  DH*DT 

RETURN 

END 

SUBROUTINE  INTSUM(TFEA,TFSU,TFMO,TFDRA,TDI,TDA,TDE,TDMM,TDMA, 
•>  TDLAN.TDH.TDAP) 
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*  CALCULATE  THE  NEW  POSITION  AND  VELOCITY  VECTORS 


SUBROUTINE  NPOS(RI,RJ,RK,R,LAN,AP.INC,  TA,A,E) 

THIS  SUBROUTINE  CALCUUTES  THE  NEW  POSITION  VECTOR 

DOUBLE  PRECISION  XV.YV.ZV, INC.RI.RJ.RK.R.LAN.AP.TA.A.E 

CALCULATE  POSITION  VECTOR  IN  'PQV'  SYSTEM 
X  «  (A*(l  -  (E**2)))/( 1  +  E*DCOS(TA)) 

XV  »  R*DCOSf,TA) 

YW  -  R*DSIN(TA) 

2W  ■  0 

TRANSFORM  POSITION  TO  'IJK'  SYSTEM 
CALL  PQVfIJK(  LAN , AP , INC , XV ,  IV , 2V , RI , R J , RK) 

R  «  DSQRT((RI**2)  +  (RJ**2)  +  (RK**2)) 

RETURN 

END 


SUBROUTINE  NVELCE.P.TA.LAN.AP.INC.VI.VJ.VK.V.HU) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  VELOCITY  VECTOR 

DOUBLE  PRECISION  INC.VP.VQ.W.MU.E.P.TA.LAN.AP.VI ,VJ,VK\V 

CALCULATE  VELOCITY  IN  'PQV'  SYSTEM 
VP  *  DSQRT(MU/P)*( -DSIN(TA) ) 

VQ  «  DSQRT(MU/P)*(E  +  DCOS(TA)) 

W  «  0.  OD+OO 

TRANSFORM  VELOCITY  INTO  ' IJK'  SYSTEM 
CALL  FQVI JK( LAN , AP , INC , VP , VQ , VV , VI , VJ , VK) 

V  «  DSQRT((VI**2)  +  (VJ**2)  +(VK**2)) 

RETURN- 

END 


*  VELOCITY  CHANGE 


SUBROUTINE  CHGVEL( DT , PER , AL , LAN , AP , I , RI , R 3 , RK , R , 

+  VI  ,VJ,VK,V,MU,PI  ,H,A,E,N,TA,P,MM,MA.,tA, 

+  TF ,T, NUM , RIRAY ,RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY .TIMRAY , 

+  TT.EI.EJ.EK.LPjHI.HJ.IOPTl.TFEA.TFSU.TFMOJFDRA, 

+  TDI ,TDA ,TDE , TDMM , TDMA , TDLAN , TDK , TD A  P ) 

*  THIS  SUBROUTINE  CALCULATE  VELOCITY  CHANGES 

*  THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

*  TACHG  *  RETURNS  TRUE  ANOMALY  FOR  VELOCITY  CHANGE  LOCATION  (CHTA) 

*  AND  AN  INDICATOR  OF  LOCATION  (ITA) 

*  CALCEl  =  CALCULATE  Orbital  ELEMENTS 

*  UN PEE!  =  CALCULATE  UNPERTURBED  ORBIT 
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*  NFQS  *  CALCULATE  NEW  POSITION 

*  NVEL  «  CALCULATE  NEW  VELOCITY 

*  STORE  *  STORE  POSITION  AND  ELEMENTS  IN  ARRAYS 

*  ENERGY  *  ENERGY  OF  SATELLITE 

*  ECC  *  ECCENTRICITY 

*  SMAXTS  «  SEMI -MAJOR  AXIS 

DOUBLE  PRECISION  T.DT.PER.AL.LAN.AP.I.RI.RJ^K.R.VI.VJ.VK.V, 

+  MU,PI ,H,A,E,N,TA,P,MM,MA,EA,TF,TT, 

+  NEWVI ,  NEWVJ ,  NEWVK ,  NEWV ,  VMAX ,  C1STA ,  E I ,  E  J ,  EK ,  LP ,  K I ,  H  J ,  VCI R , 

+  DI , DE , DA , DMM , DMA , DLAN , DH , DAP » NAVE I , NEVE J , NEWER , NEVE , NEVENR , 

+  NEWA,NEVR?,RE 

DIMENSION  RARAY(SOO) ,TARAY(500) .RIRAY(SOO) ,RJRAY(500) , 

+  RKRAY(SOO) ,AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTERS 1 , YORN , PYORN 

RE  »  6. 3782D+03 

*  PROMPT  THE  USER  FOR  TH,;  '’ELOCITY  Change  LOCATION 
CALL  TACNG( PI ,CHTA , ITA) 

*  SET  TIME  COUNTER  TO  ONE  TIME  STEP 
T  =  DT 

*  ROTATE  TO  THE  VELOCITY  CHANGE  LOCATION 

*  THIS  IS  IDENTICAL  TO  THE  Unperturbed  ORBIT  WITH  THE  EXCEPTION 

*  THAT  A  COMPLETE  ORBIT  IE  NOT  CALCULATED 
PRINT*,’ ROTATE  TO  VELOCITY  CHANGE  LOCATION’ 

IF  ((ITA.  EQ.  2)  .OR.  (ITA.EQ.3))  THEN 

PRINT*, ’BEFORE  TA  »’ ,TA 
IF  (TA  .  GT.  6.21)  THEN 
TA  «  TA  -  (2*PI) 

END  IF 

250  IF((T.  LE.  PER).  AND.  (TA.  LT.  CJiTA))  THEN 

*  PRINT*, ’TA  *’ ,TA 
NUN  =  NUN  +  1 

TT  *  TT  +  DT 

CALL  NEWELT(MM , MA , E , EA ,TA ,TF , DT, PI , PER) 

CALL  NPOS(RI ,RJ,RK,R,LAN,AP, I ,  TA,A,E) 

CALL  NVEL(E,P,TA,LAN,AP, I ,VI ,VJ,VK,V,MU) 

CALL  STORE(RI , R J , RK , R , TA , RIRAY , RJRAY , RKRAY , RARAY , 

+  TARAY , NUM , I , AP , A I NRAY , APRA Y , TT , TIMRAY ) 

T  =  T  +  DT 
GOTO  250 

ENDIF 

IF  (TF  .GE.  PER)  THEN 
TF  =  TF  -  PER 

ENDIF 

ENDIF 

*  PRINT  ESCAPE  VELOCITY  AND  CIRCULAR  VELOCITY  FOR  Reference 
CALL  EXCMS(’CLRSCRN’) 

PRINT*, ’AFTER  TA  =’ ,TA 

PRINT*, ’THIS  SHOULD  BE  THE  DESIRED  RADIUS  RP  OR  RA’ 
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260  PRINT*, 'RADIUS  ■' ,R 
PRINT'S 'VELOCITY  *'  ,V 
VMAX  =  DSQRT( 2. 0*(MU  /  R)) 

PRINT'S*. MAX  VELOCITY  AT  THIS  RADIUS  IS:1 , VMAX 
VC IS  =  DSQRTCMU/R) 

PRINTS t 'CIRCULAR  VELOCITY  AT  THIS  RADIUS  IS  : ' ,VCIR 

*  PROMPT  USER  TO  CHANGE  VELOCITY  IN  ORBITAL  PLANE 

PRINT*, ‘DO  YOU  WANT  TO  CHANGE  THE  VELOCITY  IN  THE  ORBITAL  PLANE?' 
PRINT*, 'ENTER  "Y"  OR  "N"  :  ' 

READ* , PYORN 

PRINT*, PYORN 

IF  (PYORN  .EQ.  'Y')  THEN 

PRINT*, 'GIVE  THE  TOTAL  CHANGE  IN  VELOCITY,  I.E.  5.0  KM.' 
PRINT*, 'THE  PROGRAM  WILL  FIGURE  OUT  THE  FINAL  VELOCITY  VECTOR' 
PRINT'S'  ENTER  VELOCITY  CHANGE:  ' 

READ*, CHGV 
PRINT'S  CHGV 

*  CALCULATE  NEW  VISOCITY  FOR  CHANGE  IN  THE  ORBITAL  PLANE 

NEWVI  *  VI  +  (CHGV  *  VI  /  V) 

NEWVJ  "■  VJ  +  (CHGV  *  VJ  /  V) 

NEWVK  «  VK  +  (CHGV  *  VK  /  V) 

*  Velocity  CHANGE  OUT  OF  ORBITAL  PLANE 
ELSEIF  (PYORN  .  EQ.  'N')  THEN 

PRINT'S1  ENTER  THE  NEW  VELOCITY  VECTOR:  ' 

PRINTS'  ENTER  THE  NEW  VI' 

READ*, NEWVI 
PRINT'S  NEWVI 

PRINTS'  ENTER  THE  NEW  VJ' 

READ*, NEWVJ 
PRINT'S  NEWVJ 

PRINT,'  ENTER  THE  NEW  VK' 

READ*, NEWVK 
PRINT"' , NEWVK 
NUM  =  1 
ITA  »  3 
ELSE 

CALL  EXCMS( ' GLRSCRN ' ) 

GOTO  260 
ENDIF 

*  PRINT  NEW  VELOCITY  FOR  USER  TO  CHECK 

NEWV  *  DSQRT( ( NEWVI**2 )  +  (NEWVJ**2)  +  (NEWVK**2)) 

PRINT, 'NEW  VI  =' ,  NEWV  I 
PRINT, 'NEW  VJ  =' , NEWVJ 
PRINT, 'NEW  VK  , NEWVK 
PRINT,' NEWV  ='  ,NEWV 

PRINT,'  ARE  THESE  VALUES  THE  ONES  YOU  WANT?' 


OR 


"N" 


PRINT'S 'ENTER  "Y" 

READ*, YORK 
PRINT,  YORN 

IF  (YORN  .EQ.  'N')  THEN 
CALL  EXCMS('CLRSCRN') 
GOTO  260 
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END  IF 

*  CHECK  FOR  VALID  VELOCITY 
IF  (  NEW  .GT.  VMAX)  THEN 

PRINT*,' YOUR  VELOCITY  IS  TO  GREAT  !!  1 
GOTO  260 
ENDIF 

*  Calculate  PERIGEE  RADIUS  TO  SEE  IF  SATELLITE  WILL  IMPACT  EARTH 
CALL  ENERGYC NEHV , R , MU , NEVENR) 

CALL  ECC( R I , R J , RK , R , NEW V I , NEW V J , NEWVK , NEWV , NEVE I , NEVE J, NEVER, 

+  NEVE, MU) 

CALL  SMAXISCMU, NEVENR, NEVA) 

NEVRP  «  NEVA*( 1. 0  -  NEVE) 

IF  (NEVRP  .LE.  RE)  THEN 

PRINT*,' YOUR  VELOCITY  AT  THIS  POINT  IS  TO  SMALL!!!  ' 
PRINT-v.'THE  SATELLITE  VILL  IMPACT  THE  EARTH!!  ' 

PRINT'1', 'THE  SATELLITES  RADIUS  OF  PERIGEE  VOULD  BE  ', NEVRP 
PRINT*,' A  NEV  VELOCITY  VILL  HAVE  TO  BE  ENTERED!!  ’ 

GOTO  260 
ENDIF 

*  ACCEPT  NEV  VELOCITY 
VI  *  NEWT 

VJ  *  NEVVJ 
VK  *  NEWS 
V  *  NEW 

*  CALCULATE  NEV  ELEMENT  VITH  NEV  VELOCITY  AND  SET  TIME  STEP 
CALL  CALCEL(  RI ,RJ,RK,R,VT ,VJ,VK,V,EI ,EJ,EK,E,A,I ,I»AN,LP,TA, 

+  PER,EA,MA,AP,ALfTF,P,PI,MU,MM,N,H,HI,HJ) 

DT  *  PER/50. 0 
T  *  DT 

*  THE  FOUR  Different  CASES  OF  VELOCITY  CHANGES  FOLLOVS: 

*  VELOCITY  CHANGE  AT  PERIGEE,  AND  NEW  >  V  Circular 

IF( ( ITA.  EQ.  1).  AND.  (NEW.  GT.  VCIR) )THEN 

CALL  UNPRET(DT,PER,AL,LAN,AP,I ,RI ,RJ,RK,R,VI ,VJ,VK,V, 

+  MU,PI,H,A,E,N,TA,P,.MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY,RKRAY,RARAY, 

+  TARAY , AINRAY , APRAY ,TIMRAY ,TT) 

*  Change  VELOCITY  AT  PERIGEE,  AND  NEW  <«  V  CIRCULAR 

*  APOGEE  AND  PERIGEE  SVAP 

ELSEIF  ( ( ITA.  EQ.  1).  AND.  (NEW.  LE.  VCIR) )THEN 

*  CLEAR  PREVIOUS  PLOTS 
NUM  =  1 

CALL  STORE( RI , R J, RK , R ,TA , RIKAY , RJRAY , RKRAY , RARAY , TARAY, 

•f  NUM ,  I ,  AP ,  AINRAY ,  APRAY ,  IT ,  TIMRAY) 

T  =  PER/2  , 

*  STEP  SATELLITE  TO  NEV  PERIGEE,  ONLY  A  HALF  ORBIT 
CALL  UNPRET( DT , PER , AL , LAN , AP , I , RI , R J , RK , R , VI , VJ , VK , V , 

+  MU, PI ,H,A,E,N,TA,P,MM, 
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+  MA,EA,TF,T,NUM,RIRAY,RJRAY, RKRAY, RARAY, 

t  TARAY , AINRAY , APRAY .TIMRAY ,TT) 

•r  RESET  TIME  COUNTER  TO  ONE  TIME  STEP 

T  *  DT 

CALCULATE  COMPLETE  NEXT  ORBIT 

CALL  UNPRETC DT,PER,AL,LAN,AP,I,RI ,RJ,RK,R,VI ,VJ,VK,V, 

+  MU.PI ,H,A,E,N,TA,P,MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY,RKRAY, RARAY, 

+  TARAY, AINRAY, APRAY, TIMRAY.TT) 

*  CHANGE  VELOCITY  AT  APOGEE,  AND  NEW  V  <  V  CIRCULAR 
ELSEIF  ((ITA.EQ.2)  .  AND.  (NEWV  .  LT.  VCIR))  THEN 

T  «  PER/2 

*  FINISH  ORBIT 

CALL  UNPRET( DT, PER, AL, LAN, AP, I, RI ,RJ ,RK,R,VI,VJ,VK,V, 

+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY, RKRAY, RARAY, 

+  TARAY , AINRAY , APRAY ,TIMRAY ,TT) 

*  CHANGE  VELOCITY  AT  Apogaa,  AND  NEWV  >■  V  CIRCULAR 

*  OR  AT  ANY  OTHER  TRUE  Anomaly 

ELSEIF  ((( ITA.EQ.2).  AND.  (NEWV.  GE.  VCIR))  .OR.  (ITA.EQ.3))  THEN 

IF  (TA  .GT.  6.21)  THEN 
TA  «  TA  -  (2*PI) 

ENDIF 

*  CLEAR  PREVIOUS  ORBITS  AND  STEP  SATELLITE  TO  NEW  PERIGEE 

T  *  TF 

NUN  =  1 

CALL  STORE( RI , R J ,  RK ,  R  ,TA ,  KIP.AY ,  RjRAY , RKRAY ,  RARAY , TARAY , 

+  NUN , I , AP .AINRAY , APRAY .TT.TIMRAY) 

CALL  UNPRETC  DT,  PER,  AL,  LAN,  AP,  I  ,P.I  ,RJ,RK,R,VI,VJ,VK,V, 

+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY, RKRAY, RARAY, 

+  TARAY, AINRAY, APRAY, TIMRAY.TT) 

IF  (TF  .GE.  PER)  THEN 
TF  *  TF  -  PER 

ENDIF 

*  CALCULATE  COMPLETE  NEXT  ORBIT 

T  *  DT 

CALL  UNPRET(DT,PER,AL,LAN,AP, I ,RI ,RJ ,RK,R,VI,VJ,VK,V, 

+  MU,PI ,H,A,E,N,TA,P,MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY, RKRAY, RARAY, 

+  TARAY, AINRAY, APRAY, TIMRAY.TT) 

ENDIF 

RETURN 

END 

*»WrVttfnViV»1tiV»ViWtiVVnViV>VVc»ViVirVr,VVnV*A*^V>Wf*,Wt*^r*,VVh>»ViV»Wlr*AjV*»V******MViV*»V.n'nV 


0RB22330 

0RB22340 

0RB22350 

0RB22360 

0RB22370 

0RB223B0 

0RB22390 

0RB22400 

0RB22410 

ORB22420 

0RB22430 

0RB22440 

0RB224S0 

ORB22460 

0RB22470 

0RB224S0 

ORB22490 

ORB22500 

0RB22S10 

ORB22520 

ORB22530 

0RB22540 

ORB22SSO 

ORB22S60 

ORB2257C 

0RB22580 

0RB22590 

ORB22600 

0RB22610 

0RB22620 

0RB22630 

ORB22640 

0RB22650 

0RB22660 

0RB22670 

ORB22680 

0RB22690 

ORB22700 

0RB22710 

0RB22720 

ORB22730 

ORB22740 

0RB22750 

0RB22760 

0RB22770 

0RB22780 

0RB22790 

ORB22800 

ORB22810 

ORB22820 

0RB22830 

ORB22840 

0RB22850 

0RB22860 

0RB22870 

0RB22880 
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SUBROUTINE  TACNGCPI , CHTA, ITA 1 

THIS  SUBROUTINE  Asks  THE  USER  FOR  VELOCITY  CHANGE  LOCATION 
DOUBLE  PRECISION  CHTA, PI 
CALL  EXCMS('CLRSCRN') 

PRINT'S 'WHERE  DO  YOU  WANT  TO  CHANGE  THE  VELOCITY?' 

PRINT'S’  I.  AT  CURRENT  PERIGEE' 

PRINT*,'  2.  AT  CURRENT  Apogat' 

PRINT'S'  3.  AT  A  SPECIFIC  TRUE  Anomaly' 

PRINT'S 'ENTER  "l",  "2"  OR  "3"' 

READ’S  ITA 
PRINTS  ITA 

SET  TRUE  ANOMALY  CHANGE  LOCATION  (CHTA)  TO  DESIRED  LOCATION 
IF  (ITA  .EQ.  1)  THEN 
CHTA  *0.0 
ENDIF 

IF  (ITA  .EQ.  2)  THEN 
CHTA  *  PI 
ENDIF 

IF  (ITA  .EQ.  3)  THEN 

PRINTS 'AT  WHAT  TRUE  ANOMALY  DO  YOU  WANT  TO  CHANGE  THE' 
PRINT'S 'VELOCITY?' 

PRINT’S 'ENTER  TRUE  ANOMALY  IN  DEGREES' 

READ’S  CHTA 
PRINTS  CHTA 
CHTA  *  CHTA  *  PI  /  ISO 
ENDIF 
RETURN 
END 


*  OUTPUT  FLOTS 

VruVnV.V 


SUBROUTINE  PLOTS( R IRAY , RJRAY , RKRAY , RARAY ,TARAY , NUM , PI , INC , LP , A , 

+  E,TF, AINRAY , APRAY ,TIMRAY ,TFEA ,TFSU ,TFMO ,TFDRA , 

+  PER ,TDI . TDA ,TDE , TDMM , TDM A , TDLAN ,TDH,TDAP , 

+  MM,MA,LAN,H,AP,R,V) 

*  THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  TYPE  OF  OUTPUT  THAT  IS 

*  DESIRED  PERIFOCAL,  GROUND  TRACK  OR  TO  SKIP  THE  PLOT. 

*  THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

*  PERIF  =  PLOT  PERIFOCAL  ORBIT 

*  GRTRK  =  PLOT  GROUND  TRACK 

*  DATE  =  DISPLAYS  DATA  ON  PLOT 

*  TEC618  *  SET  Disspla  TO  TEC  618  OUTPUT 

*  ENDPL  *  END  THIS  DISSPLA  PLOT 

*  REFER  TO  DISSPLA  USER"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 

*  SUBROUTINES 

DOUBLE  PRECISION  PI , A , E , INC , LP ,TF , PER , MM , MA , LAN , H , AP , R , V 
DIMENSION  RIRAY(500) ,RJRAY(500) ,RKRAY(500) ,RARAY(500) ,TARAY(500) , 


0RB22890 

ORB22900 

QRB22910 

0RB22920 

0RB22930 

0RB22940 

0RB22950 

0RB22960 

0RB22970 

ORB22980 

ORB22990 

0RB23000 

ORB23010 

ORB2302O 

ORB23030 

ORB23040 

ORB23050 

ORB23060 

QRB23070 

ORB23080 

ORB23090 

0RB23100 

0RB23110 

0RB23120 

0RB23130 

0RB23140 

ORB23150 

0RB23160 

ORB23170 

0RB23I80 

0RB23190 

ORB23200 

ORB232IO 

0RB23220 

ORB23230 

0RB23240 

0RB23250 

ORB23260 

0RB23270 

0RB23280 

0RB23290 

ORB23300 

ORB23310 

0RB23320 

ORB23330 

ORB23340 

ORB23350 

ORB23360 

0RB23370 

ORB23380 

ORB03390 

ORB23400 

0R323410 

ORB23420 

ORB23430 

ORB23440 
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+  AIN8AY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* 1 ,YORN 

CALL  EXCMSC 1 CLRSCRS' ) 

*  CALCULATE  SINGLE  PRECISION  VARIABLES 

SPI  «  SNGL(PI) 

SA  *  SNGL(A) 

SE  «  SNGL(E) 

SINC  -  SNGL(INC) 

SLP  *  SNGL(LP) 

STF  *  SNGL(TF) 

SPER  «  SNGL(PER) 

SMM  *  SNGL(MM) 

SMA  -  SNGL(MA) 

SLAN  -  SNGL(LAN) 

SH  *  SNGL(H) 

SAP  «  SNGL(AP) 

SV  «  SNGL(V) 

SR  *  SNGL(R) 

*  PROMPT  USER  FOR  DISPLAY  TYPE 

340  PRINT*, 'WHAT  TYPE  OF  DispUy  IS  DESIRED: ' 

PRINT*, '  1.  PERIFOCAL’ 

PRINT*, '  2.  GROUND  TRACK' 

PRINT*, '  3.  SKIP  PLOT' 

PRINT*, 'ENTER  1,2, 3, 4: ' 

READ*, INPUT 
PRINT350, INPUT 
350  FORMAT(IA) 

CALL  TEK618 

*  CALL  APPROPRIATE  PLOT 
IF  (INPUT  .EQ.  1)  THEN 

CALL  PERI F( RARAY .TARAY , NUM , SPI , S INC , SLP , SA , SE) 

ELSEIF  (INPUT  .  EQ.  2)  THEN 

CALL  GRTRK( AINRAY , APRAY , TARAY , STF , NUM , TIMRAY) 

ELSE7F  (INPUT  .  EQ.  3)  THEN 
GGfO  360 
ELSE 

PRINT*, ' INVALID  ENTRY! ' 

GOTO  340 
FNDIF 

*  DISPLAY  DATA 

CALL  DATA( SINC , SA , SE ,TFEA ,TFSU ,TFMO ,TFDRA , SPER , SPI ,TDI ,TDA ,TDE , 
+  TDMM ,  TDMA ,  TDLAN ,  TDH ,  TDAP ,  SMM ,  SMA ,  SLAN ,  SH ,  SAP ,  SV ,  SR  ) 

CALL  ENDPL(O) 

*  PROMPT  USER  IF  ANOTHER  DISPLAY  TYPE  IS  DESIRED 

PRINT*, 'WOULD  YOU  LIKE  ANOTHER  PLOT  USING  THE  SAME  ORBITAL' 
PRINT* , ' PARAMETERS  AND  DATA: ' 

PRINT* , ' ENTER  "Y"  OR  "N"  :' 

READ* , YORN 
PRINT*, YORK 


0RB23450 

0RB23460 

0RB23470 

0RB23480 

0RB23490 

0RB23500 

0RB23510 

ORB23520 

0RB23530 

0RB23540 

ORB23550 

0RB23560 

ORB23570 

ORB235BO 

ORB23590 

ORB23600 

ORB23610 

0RB23620 

0RB23630 

0RB23640 

ORB23650 

0RB23660 

0RB23670 

0RB236B0 

0RB23690 

0RB23700 

0RB23710 

0RB23720 

0RB23730 

ORB23740 

0RB23750 

0RB23760 

0RB23770 

0RB23780 

0RB23790 

0RB23800 

ORB23810 

0RB23820 

ORB23830 

ORB23840 

0RB23850 

0RB23660 

0RB23870 

0RB23880 

0RB23890 

ORB23900 

0RB23910 

0RB23920 

0RB23930 

ORB23940 

0RB23950 

0RB23960 

ORB23970 

0RB23980 

0RB23990 

ORB24000 


6$ 


IF  (YORK  .EQ.  *Y')  THEN’ 
GOTO  340 
ESDIF 

360  RETURN 
END 


SUBROUTINE  PERIF( RARAY ,TARAY , NUM , PI , INC , LP , A , E ) 

*  THIS  SUBROUTINE  PLOTS  OUT  THE  RESULTS  OF  THE  PROGRAM  USING  THE 

*  DISPLAY  FEATURE  ON  THE  MAIN  FRAME. 

*  REFER  TO  DISSPLA  USERS  GUIDE  FOR  EXPLANATION  OF  DISSPLA 

*  SUBROUTINES. 

REAL  INC,LP 

DIMENSION  TARAY(500) ,RARAY(500) ,RIRAY(500) ,RJRAY(500) ,RKRAY(500) 
I  *  1 

*  SET  SCALE  OF  AXIS 
RSTEP  «  (A,v(  1+E))  /  3 
CALL  TEN61S 

CALL  RESET( 3HALL) 

CALL  SCMPIX 

CALL  PHYSORC  1.25,4.  ) 

CALL  AREA2D(6. ,6. ) 

CALL  MESSAGC ’ PERIFOCAL  COORDINATE  SYSTEM? 1 ,100,1.0,6.5) 

CALL  XNANS<'XV\2) 

CALL  YNAME('YV',2) 

CALL  XAXANG(90. 0) 

CALL  YAXANG(O.O) 

CALL  INTAXS 

CALL  POLAR( 1.  , RSTEP, 3.  ,3.  ) 

CALL  POLY 3 
CALL  NOCHEK 

CALL  CURVE(TARAY, RARAY, NUM, 1) 

CALL  COMPLY 
CALL  HSIGHTC.2) 

CALL  RSSET(' COMPLEX') 

CALL  RESETC HEIGHT') 

CALL  ENDGR(O) 

*  Display  EARTH  PLOT 

CALL  EARTH1(A,E, INC, PI, RSTEP) 

RETURN 

END 

IrSt'lrlt'Mt'kititMtMffrMrMtltMrMtirMtMtitMtititMtMtMrMtltitiitlrMtirtsirMt 

SUBROUTINE  EARTHlC A , E , INC , PI , RSTEP) 

*  THIS  SUBROUTINE  PLOTS  A  VIEW  OF  THE  WORLD,  LOOKING  DOWN  THE  'Z' 

*  AXIS,  PLACED  ON  THE  ORIGIN.  THE  Latitude  IS  FIXED,  BUT  THE 

*  LONGITUDE  VARIES  WITH  THE  INCLINATION. 

*  REFER  TO  DISSPLA  USER"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 

*  SUBROUTINES 


0RB24010 

0RB24020 

ORB24030 

ORB24040 

0RB24050 

ORB24060 

0KB24070 

ORB240B0 

ORB24090 

0RB24100 

0RB24U0 

URB24120 

ORB24130 

0RB24140 

ORB24150 

ORB24160 

ORB24170 

QRB24180 

ORB24190 

0RB24200 

0RB24210 

0RB24220 

0RB24230 

0RB24240 

0RB24250 

0RB24260 

0RB24270 

0RB24280 

ORB24290 

ORB24300 

ORB24310 

ORB24320 

ORB24330 

0RB24340 

0RB24350 

ORB24360 

0RB24370 

ORB24380 

ORB24390 

ORB244CO 

0RB24410 

0RB24420 

0RB24430 

0RB24440 

0RB24450 

0RB24460 

ORB24470 

ORB24480 

0RB24490 

ORB24500 

0RB24510 

0RB24520 

ORB24530 

0RB24540 

ORB24550 

0RB24560 
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REAL  INC 

COMMON  IV0RKO8C0) 

DATA  IVDIM/3500/ 

RE  ■  6375. 145 

*  SCALE  THE  EARTH  PLOT  AND  CENTER  ON  THE  ORIGIN 
SCFAC  ■  RE/RSTEP 

SCFAC2  -  SCFAC  *  2.  0 
XPHS  *  1.  25  +  3.  0  -  SCFAC 
YPHS  »  4.  0  +  3.  0  -  SCFAC 
YPOLE  »  90  -  (INC  *  180  /  PI) 

IF(YPOLE  .GT.  90)  THEN 
YPOLE  «  YPOLE  -  90 
ENDIF 

YORIG  «  YPOLE  -  90 
YMAX  »  YPOLE  +  90 
CALL  RESETOHALL) 

CALL  PHYSOR( XPHS. YPHS) 

CALL  PROJCTC LAMBERT  EQ/AREA') 

CALL  MAPOLEC 0.0 .YPOLE) 

CALL  AREA2D  (SCFAC2.SCFAC2) 

CALL  THKFRM  (0. 02) 

CALL  GRAF(-90.  ,30.  ,90.  , YORIG, 30.  ,YMAX) 

CALL  FRAME 

CALL  MAPFIL(’HAPDTA') 

CALL  LBLANKC 'LAND1 ,IWDIM) 

CALL  GRIDC1, I> 

CALL  LBLANKC 'WATER' , IVDIM) 

CALL  DASH 
CALL  GRID(1,1) 

CALL  RESET  ('DASH') 

CALL  *,NDGR(0) 

RETUR. 

END 

SUBROUTINE  GRTRK(AINRAY.,APRAY,TARAY,YF,NUM,TIMRAY) 

DIMENSION  AINRAYC 500)  ,AP1'AY(500)  ,TARAY(500)  , 

+  ELARAY(500) ,ELORAY(500) ,TLONG(500) ,TLAT(500) ,TIMRAY(500) 

RE  »  6. 3782E+03 
EROT  »  7. 292115856E-05 
STF  *  (TF) 

I  «  1 

*  LOAD  ARRAYS  WITH  LATITUDE  AND  LONGITUDE 
410  IF  (I  .LE.  NUM)  THEN 

X  =  RE*CQS(APRAY(I))*COS(TARAY(I))-RE*SIN(APRAY(I))* 

+  SIN(TARAY(I)) 

Y  =  RE*COS( AINRAYC I ) )*SIN( APRAYC I ) )*C0S(TARAY( I) )  + 

+  R2*C0S( AINRAYC I ) )*COS( APRAY( I ) )*SIN(TARAY( I ) ) 

Z  =  RE*SIN( AINRAYC I) )*SIN( APRAYC I) )*COS(TARAY( I) )  + 

+  RE*SIN( AINRAYC I ) )*COS( APRAYC I) )*SIN(TARAY( I ) ) 


0RB24570 
0RB24580 
0RB24590 
0RB24600 
0RB24610 
0RB24620 
ORB 246 30 
0RB24640 
0RB24650 
0RB24660 
0RB24670 
ORB24680 
0RB24690 
0RB24700 
0RB24710 
ORB24720 
ORB24730 
ORB24740 
0RB24750 
0RB24760 
0RB24770 
0RB24760 
0RB24790 
0RB24800 
0RB24810 
0RB24820 
0RB24830 
0RB24840 
ORB24850 
0RB24860 
0RB24870 
0RB24880 
ORB24890 
ORB24900 
0RB24910 
ORB24920 
ORB24930 
ORB24940 
ORB24950 
0RB24960 
0RB24970 
ORB24980 
0RB24990 
ORB25000 
0RB25010 
ORB25020 
0RB25030 
0RB25040 
ORB25050 
0RB25060 
ORB25070 
ORB25080 
0RB25090 
ORB25100 
0RB25110 
0RB25120 
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*  CALCULATE  LATITUDE 

ELARAY( I )  *  (ASIX(Z/RE))  *  (160/3.14155) 

*  TEA?  'X'  AND  'Y'  FOR  ARCTAN  IX  CALCULATING  LOXGITUDE 
IFUY  .LE.  10)  .AND.  (Y  .GE.  0.0))  THEN* 

Y  «  10. 

ELSEIF  ((Y  .  GE.  *10).  AND.  (Y  .  LE.  0.0))  THEN 

Y  «  -10. 

ENDIF 

IF((X  .LE.  10)  .  AXD.  (X  .GE.  0.0))  THEN 
X  *  10. 

ELSEIF  ((X  . GE.  -10).  AND. (X  .LE.  0.0))  THEN 
X  -  -10. 

ENDIF 

*  CALCULATE  LONGITUDE 

ELORAY(I)  *  (ATAN2(Y,X)  -  ( EROT*TIMRAY( I ) ) )  *  (180/3.14159) 

*  MODIFY  LONGITUDES  TO  (  -ISO  TO  ISO) 

420  IF  (ELORAY(I)  .  LT.  -ISO)  THEN 

ELORAY(I)  «  ELORAY(I)  +  360 
GOTO  420 

ENDIF 
I  -  I  +  1 
GOTO  410 
ENDIF 

*  SET  DISSPLA 
CALL  TSK61S 
CALL  RESET( SHALL) 

CALL  YAXANG  (0. ) 

CALL  PHYS0R(1. 0,6.0) 

CALL  XXAMEC'  ',1) 

CALL  YNAME(‘  ',1) 

CALL  AREA2D(7.  5 ,3.  75) 

CALL  HEADIN  ('GROUND  TRACKS' ,100,1- 5,1) 

CALL  SCMPLX 

CALL  MA?GR( -130.  ,90.  ,180.  ,-90.  ,30.  ,90.) 

CALL  GRID  (1,1) 

CALL  MAPFIL  ( 'MAPDTA* ) 

I  «  1 

*  IGNORE  Boundary  POINTS 

430  IF  ((ELORAY(I)  .  LT.  -175)  .OR. 

+  (ELORAY(I)  .GT.  175)  .OR. 

+  (ELARAY(I)  .LT.  -85)  .OR. 

+  (ELARAY(I)  .GT.  85))  THEN 

1  =  1  +  1 
GOTO  430 
ENDIF 


ITEMP  =  1 

LOAD  FIRST  POINT  OF  NEW  PLOT  SEGMENT 
IF  (I  .LE.  NUM)  THEN 


0RB25130 

0RB25140 

QRB2S150 

0RB25160 

0RB25170 

0RB25180 

0RB25190 

0RB25200 

0RB25210 

0RB25220 

0RB25230 

0RB25240 

0RB25250 

0RB25260 

0RB25270 

0RB2S280 

0RB25290 

ORB25300 

ORB25310 

0RB25320 

0RB25330 

0RB25340 

0RB25350 

0RB25360 

0RB2S370 

0RB25380 

0RB25390 

0RB25400 

0RB2541Q 

0RB25420 

0RB25430 

0RB25440 

0RB25450 

0RB25460 

0RB25470 

0RB25480 

ORB25490 

0RB25500 

0RB25510 

ORB25520 

0RB25520 

0RB25540 

ORB25550 

0RB25560 

ORB25570 

0RB25580 

0RB25590 

ORB25600 

ORB25610 

ORB25620 

ORB25630 

ORB25640 

ORB25650 

ORB25660 

0RB25670 

ORB25680 


TLONG( ITEM?)  «  ELORAY(I) 

ORB25690 

TLAT( ITEM?)  «  ELARAV(I) 

0RB2S700 

1  =  1+1 

0RB2S710 

Vr 

IF  (  I  .GE.  SUM)  THEN 

0RB2S720 

* 

CALL  F0LV3 

0P.B25730 

tV 

CALL  CURVE( TLOSG ,TLAT , ITEMP , 1 ) 

0RB2S740 

* 

ENDIF 

ORB2S7SO 

END  IF 

0RB2S760 

ORB2S770 

* 

LOAD  SECOND  POINT  IN  LINE  SEGMENT 

ORB25780 

IF  (I  .LE.  Nl’M)  THEN 

0RB25790 

ITEMP  «  ITEMP  +  1 

QRB2S800 

TLONG(ITEMP)  -  ELORAY(I) 

0RB2S810 

TLAT( ITEMP)  -  ELARAY(I) 

0RB2S820 

1  =  1  +  1 

0RB25830 

IF  (  I  .GE.  NUM)  THEN 

ORB 25 840 

CALL  P0LY3 

0RB25850 

CALL  NOCHEK 

0RB25860 

CALL  CURVE(TLONG,TLAT, ITEMP, 1) 

0RB2S870 

ENDIF 

0RB25880 

ENDIF 

0RB2S890 

0RB25900 

★ 

LOOP  UNTIL  SEGMENT  REACHES  EDGE  OR  NO  MORE  POINTS 

0RB2S910 

440 

IF  (I  .LE.  NUM)  THEN 

0RB25920 

0RB25930 

* 

BOTH  LAT  AND  LONG  INCREASING 

0R325940 

IF((ELORAY(I  -  2)  .  LE.  ELORAYCI  -  1))  .  AND. 

0RB25950 

+  (ELARAYC I  -  2)  .  LE.  ELARAYCI  -  1)))  THEN 

0RB25960 

IF((ELGRAY(I)  .  LT.  -170)  .OR. 

0RB25970 

+  (ELARAYCI)  .LT.  -80))  THEN 

0RB25980 

CALL  POLY 3 

0RB25990 

CALL  NOCHEK 

ORB26000 

CALL  CURVE(TL0NG ,TLAT, ITEMP , 1) 

0RB26010 

GOTO  430 

ORB26020 

ELSE 

0RB26030 

ITl  «  ITEMP  +  1 

ORB26040 

TLOMi( ITEMP)  ■  ELORAYCI) 

0RB26050 

TLAT( ITEMP)  =  ELARAY(I) 

0RB26060 

ENDIF 

0RB26070 

0RB26080 

★ 

BOTH  LAT  AND  LONG  DECREASING 

0RB26090 

ELSEIFC (ELORAYCI  -  2)  .GT.  ELORAYCI  -  1))  .AND. 

ORB26100 

+  (ELARAYCI  -  2)  .  GT.  ELARAYCI  -  1)))  THEN 

0RB26110 

IF( (ELORAYCI)  .GT.  170)  .OR. 

ORB26120 

+  (ELARAY(I)  .GT.  80))  THEN 

0RB26130 

CALL  POLY 3 

0RB26140 

CALL  NOCHEK 

0RB26150 

CALL  CURVECTLONG.TLAT, ITEMP, 1) 

0RB26160 

GOTO  430 

ORB26170 

ELSE 

0RB26180 

ITEMP  =  ITEMP  +  1 

0RB26190 

TLONGC ITEMP)  «  ELORAY(I) 

ORB26200 

TLATC ITEMP)  =  ELARAY(I) 

ORB26210 

ENDIF 

0RB26220 

0RB26230 

ii 

LAT  INCREASING,  LONG.  DECREASING 

0RB26240 

ELSE I F( ( ELOR AV( I  -  2)  .GT.  ELORAY( I  -  l))  .AND. 
+  (ELARAY(I  -  2)  .  LE.  ELARAYd  -  1)))  THEN 

IF((ELORAY(I)  .GT.  170)  .QR. 

+  l  ELARAYd)  .LT.  -80))  THEN 

CALL  P0LY3 
CALL  NOCHEK 

CALL  CURVE ( TLONG ,TLAT , ITEMP , 1) 

GOTO  430 

ELSE 

ITEMP  *  ITEMP  +  1 
TLONG( ITEMP)  -  ELORAY(I) 

TLAT( ITEMP)  *  ELARAY(I) 

ENDIF 

*  LAT.  DECREASING,  LONG.  INCREASING 

ELSEIF((ELORAY(I  -  2)  .  LE.  ELORAY(I  -  1))  .AND. 
+  (ELARAYCI  -  2)  .GT.  ELARAY(I  -  1)))  THEN 

IF((ELQRAY(I)  .LT.  -170)  .OR. 

+  (ELARAY(I)  .GT.  SO))  THEN 

CALL  P0LY3 
CALL  NOCHEK 

CALL  CURVE( TLONG ,TLAT, ITEMP ,1) 

GOTO  430 

ELSE 

ITEMP  «  ITEMP  +  1 
TLONGC ITEMP)  *  ELORAY(I) 

TLAT( ITEMP)  «  ELARAY(I) 

ENDIF 

ENDIF 

IF(  I  .EQ.  NUM)  THEN 
CALL  P0LY3 
CALL  NOCHEK 

CALL  CURVE(TLONG,TLAT, ITEMP, 1) 

ENDIF 
I  »  I  +  1 
GOTO  440 
ENDIF 


CALL  POLY 3 
CALL  NOCHEK 

CALL  CURVE(TLONG,TLAT, ITEMP, 1) 


CALL  COMPLX 
CALL  HEIGHTS  2) 

CALL  THKFRM  (0.03) 
CALL  FRAME 
CALL  RESETC COMPLX') 
CALL  RESETC 'HEIGHT') 
CALL  ENDGR  (0) 

RETURN 

END 


AVrfinViiViVKKiViVKiVitVtiViViWnn^WnWiVAiVuiWnWfMWnWiMWnlnVrtiVtVilr^ViVVhfiWnWtAAuiiAiViHti'!* 


0RB26230 
0RB26260 
0RB26270 
0RB262B0 
0RB26290 
0RB26300 
0RB26310 
0RB26320 
0RB26330 
0RB26340 
0RB26350 
0RB26360 
0RB26370 
0RB26380 
0RB26390 
0RB26400 
ORB 264 10 
0RB26420 
GRB26430 
0RB26440 
0RB26450 
0RB26460 
0RB26470 
0RB26480 
0RB26490 
0RB26500 
0RB26510 
0RB26520 
0RB26530 
0RB26540 
0RB26550 
0RB26560 
0RB26570 
0RB26580 
0RB26590 
0RB26600 
0RB26610 
0RB26620 
0RB26630 
ORB26640 
0RB26650 
0RB26660 
0RB26670 
0RB26660 
CRB26690 
ORB2670O 
0RB26710 
0RB26720 
0PB26730 
0RB26740 
ORB 26 750 
0RB26760 
0RB26770 
0RB26780 
0RB26790 
ORB26800 
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* 

it 

it 


* 

★ 


* 


* 


it 


★ 


* 


SUBROUTINE  DATA( I , A , E , TFEA , TFSU , TRIO , TFDRA , PER, PI ,TOI , TDA ,TDE , 

+  TDNM.TDNA.TDLAN.TDH.TDAP.MM^.LAS.H.AP.V.R) 

THIS  SUBROUTINE  Display*  THE  ORBITAL  DATA  FOR  BOTH  THE  PERIFOCAL 
AND  THE  GROUND  TRACK  PLOTS. 

REFER  TO  DISSPLA  US£R"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 
SUBROUTINES 

REAL 

MU  «  3. 986012E+05 

CALCULATE  THE  AVERAGE  FORCES  FROM  THE  TOTAL  MAGNITUDE  OF 

FORCE  CHANGES 

AVGFE  «  TFEA/50.  0 

AVGFS  *  TFSU  /  50.  0 

AVGFM  ■  TFMO  /  50. 0 

AVGFD  -  TFDRA  /  50. 0 

CALCULATE  ORBITAL  ELEMENTS  IN  Uiabla  UNITS 
PERM  «  PER/3600 

DI  “  I  *  (180.0/PI) 

DLAN  «  LAN  *  (180.0/PI) 

DAP  «  AP  *  (180.0/PI) 

CALCULATE  Av«r*s«  CHANGE  IN  ELEMENTS  FOR  0;'JE  PERIOD 

AVGDI  «  TDI  /  50. 0 

AVGDA  «  TDA  /  50. 0 

AVGDE  -  TDE  /  50, 0 

AVGDMM  *  *i  uM>j  /  50. 0 

AVGDMA  -  TDMA  /  50. 0 

AVGLAN  «  TDLAN  /  50,  0 

AVGDH  ■  TDH  /  50. 0 

AVGDAP  -  TDAP  /  50.  0 

CALCULATE  RADIUS'S  AND  VELOCITIES 
ENR  «  ((V**2)/2)  -  (MU/R) 

RP  -  A*(l  -  E) 

RA  «  A*( 1  +  E) 

VP  «  SQRT( 2*(ENR  +  (MU/RP))) 

VA  «  SQRT(2*(ENR  +  (MU/RA))) 


SET  DISSPLA 
CALL  RESET(3HALL) 

CALL  SCMPLX 

CALL  PHYS0R(0.  0,0.0) 

CALL  AREA2D(8.  5 ,4.  0) 

PRINT  DATA 

CALL  MESSAG('I  =  S' ,100,0. 25,3. 67) 

CALL  REALN’0(DI,3,  *ABUT' ,' ABUT') 

CALL  MESSAG( '  DEG.  $ '  100 , ' ABUT' . ' ABUT' ) 
CALL  MESSAGC  A  *  $f ,  100, '  ABUT,' ABUT' ) 
CALL  P.EALNO(A,  1,  'ABUT' , '  ABUT'  ) 

CALL  MESSAG( '  KM§ ' , 100 , ' ABUT' , ' ABUT' ) 


0RB26610 
0RB26820 
0RB26830 
0RB26840 
0RB26850 
0RB26860 
0RB26870 
0RB26880 
0RB26890 
0RB26900 
ORB 269 10 
0RB26920 
0RB26930 
0RB26940 
ORB26950 
ORB26960 
OR826970 
ORB26980 
0RB26990 
ORB27000 
0RB27010 
ORB27020 
0RB27030 
0RB27040 
0RB27050 
ORB27060 
0RB27070 
0RB27080 
0RB27090 
0RB27100 
0RB27110 
0RB27120 
0RB27130 
0RB27140 
0RB27150 
0RB27160 
ORB27170 
ORB27180 
0RB27190 
ORB 27 200 
ORB27210 
ORB27220 
ORB27230 
0RB27240 
CRB27250 
ORB27260 
0RB27270 
0RB27280 
0RB27290 
ORB27300 
0RB27310 
ORB27320 
0RB27330 
0RB27340 
0RB27350 
0RB27360 
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CALL  MESSAGC*  E  *  S' , 100, 'ABUT',' ABUT' ) 

CALL  REALNOC  E , 3 , ' ABUT' . ' ABUT' ) 

CALL  MESSAGC '  PER  *  $' ,100, 1  ABUT’ , ' ABUT' ) 

CALL  REALNOC PERM, 2, 'ABUT* , 1 A3UT* ) 

CALL  MESSAGC '  HOURS? ' , 100, ' ABUT’ , ' ABUT' ) 

CALL  MESSAGC 'AVERAGE  RATE  OF  CHANGE  OF  ELEMENTS  PER  SECOND  ?' , 

+  100,1.0,3.0) 

CALL  NHSAGC'DI/DT  -  S' ,100,0.  25,2.  67) 

CALL  REALNOC AVGDI, -2, 'ABUT' , 'ABUT' ) 

CALL  MESSAGC '  DA/DT  ■  $ ' , 100 , ' ABUT' , ' ABIT' ) 

CALL  REALNOC AVGDA , -2 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC'  DE/DT  -  $' ,100, 'ABUT' , 'ABUT' ) 

CALL  REALNOC AVGDE, -2, 'ABUT' , 'ABUT' ) 

CALL  MESSAGC 'DMM/DT  «  $', 100,0. 25,2. 33) 

CALL  REALNOC AVGDMM, -2, 'ABUT' , ' ABUT' ) 

CALL  MESSAGC '  DNA/DT  ■  5 ' , 100 , ' ABUT’ , ' ABUT' ) 

CALL  REALNOC AVGDMA , -2 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC '  DLAN/DT  -  $ ' , 100 , ' ABUT’ , ’ ABUT' ) 

CALL  REALNOC AVGLAN , -2 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC 'DH/DT  *  $' , 100,0.  25,2. 00) 

CALL  REALNOC AVGDH, -2 , ' ABUT' , ' ABUT’ ) 

CALL  MESSAGC '  DAP/DT  *  S ' , 100 , ’ ABUT' , ' ABUT' ) 

CALL  REALNOC  AVGDMA , -2 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC 'AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S**2) 
+5'  ,100,1.0,1.67) 

CALL  MESSAGC ' EARTH  *  S' ,100,0. 10,1. 33) 

CALL  REALNOC AVGFE , - 1 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC '  MOON  ■  S ' , 100 , ' ABUT' , ' ABUT' ) 

CALL  REALNOC AVGFM, -1 , 'ABUT' , 'ABUT' ) 

CALL  MESSAGC'  SUN  ■  $' ,100,  ABUT* .'ABUT') 

CALL  REALNOC AVGFS,-1,' ABUT',' ABUT') 

CALL  MESSAGC'  DRAG  -  5' ,100, 'ABUT' , 'ABUT') 

CALL  REALNOC AVGFD,-1, 'ABUT' , 'ABUT* ) 

CALL  MESSAGC ' PERIGEE? ' , 100 , 2.  73 , 1.  0) 

CALL  MESSAGC '  Apoftt? ' , 100 , ' ABUT' , ' ABUT* ) 

CALL  MESSAGC 'RADIUS  (KM)?' ,100,0. 25,0. 67) 

CALL  MESSAGC 'RP  «?' ,100,2. 75,0. 67) 

CALL  REALNOC R?,l,' ABUT', 'ABUT') 

CALL  MESSAGC'  ?', 100, 'ABUT' , 'ABUT' ) 

CALL  MESSAGC '  RA  *? ' , 100 , ' ABUT' , ' ABUT'  ) 

CALL  REALNOC RA,1,' ABUT', 'ABUT') 

CALL  MESSAGC 'VELOCITY  CKM/SEC)?'  100,0.25,0.33) 

CALL  MESSAGC ’VP  ■?' ,100,2. 75,0. V3) 

CALL  REALNOC VP, 2, 'ABUT', 'ABUT') 

CALL  MESSAGC '  S' , 100, *  ABIT' . ' ABUT* ) 

CALL  MESSAGC '  VA  *$' , 100, 'ABUT' , ‘ABUT* ) 

CALL  REALNOC  VA , 2 , ' ABUT ' , ' ABUT ' ) 
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ORB27420 

ORB27430 

0RB27440 
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ORB27540 
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0RB27600 
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ORB27630 

ORB27640 

0RB27650 

0RB27660 

0RB27670 
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0RB27690 
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0RB27720 
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0RB27750 

0RB27760 

ORB27770 

0RB27780 

0RB27790 

0RB27800 

ORB27810 

ORB27820 

0RB27830 

0RB27840 

0RB27850 

ORB27860 

0RB27870 

ORB27880 

ORB27890 

ORB27900 

0RB27910 

ORE27920 


!1WO 


CALL  RESET('COMPL>:' 
ALL  ENDGR(O) 


) 


.TURN' 


ND 


0RB27930 

0RB27940 

ORB279SO 

0RB27960 

0RB37970 

ORB279SO 
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APPENDIX  B.  COORDINATE  SYSTEMS 


A.  'UK' :  GEOCENTRIC  -  EQUATORIAL 


The  geocentric-equatorial  system  as  seen  in  Figure  3  has  its  origin  at  the  earth's 
center.  The  fundamental  plane  is  in  the  equator  and  the  positive  X-axis  points  in  the 
vernal  equinox  direction.  The  Z-axis  points  in  the  direction  of  the  north  pole.  This 
system  is  not  fixed  to  the  earth  and  turning  with  it;  rather,  the  geocentric-equatorial 
frame  is  nonrotating  with  iespect  to  th*  stars  (except  for  precession  of  the  equinoxes) 
and  the  earth  turns  relativ-  to  it.  Unit  vectors,  7,  J,  and  K  shown  in  Figure  3,  lie  along 
the  X,  Y,  and  Z  respectively,  (Ref.  1:  p.55j 
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B.  'PQW':  PERIFOCAL 


l  _  _ 

V 


Figure  4.  Perifocal  coordinate  system 


The  perifocal  coordinate  system  has  its  fundamental  plane  in  the  plane  of  the  satel¬ 
lite’s  orbit  as  seen  in  Figure  4.  The  coordinate  axis  are  named,  A'.,  and  Zm  . .  The 
A'„  axis  points  toward  the  perigee;  the  axis  is  rotated  90  degrees  in  the  direction  of 
orbital  motion  and  lies  in  the  orbital  plane;  the  Zm  axis  along  h  completes  the  right- 
handed  perifocal  system.  Unit  vectors  in  the.  direction  of  Xm,  Ym  and  Zm  are  called  P%  Q 
and  IK  respectively.  jRef.  1:  p.57) 
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'RSW' :  ORBITAL 


Figure  5.  Orbital  coordinate  system 


(Figure  9.4-1,  Ref.  1) 

The  orbital  coordinate  system  has  its  principle  axis,  R  (unit  vector  r),  along  the  in¬ 
stantaneous  radius  vector,  r  as  seen  in  Figure  5.  The  axis  S  is  rotated  90  degrees  from 
R  in  the  direction  of  increasing  true  anomaly.  The  third  axis,  W,  is  perpendicular  to 
both  R  and  S.  Note  that  this  coordinate  system  is  simply  rotated  v,  from  the  PQW 
perifocal  system.  (Ref.  1:  p.39S] 


D.  COORDINATE  TRANSFORMATIONS 

The  coordinate  transformations,  for  the  previous  coordinate  systems,  use  angular 
rotations  about  the  axis  to  evaluate  the  transformation  matrix.  The  matrix  elements  rH 
arc  calculated,  then  applied  to  the  old  vector  to  get  the  vector  in  the  new  coordinate 
system.  The  following  orbital  elements  arc  used: 

0  -  longitude  of  ascending  node 
o  *  argument  of  perigee 
i  ■  inclination 
M,  -  argument  of  latitude 
v0  ■  true  anomaly 

The  coordinate  transformations  follow  (Ref.  1:  p.74-$3) 

1.  PQWtoIJK 

r„  *  cos  £2  cos  (0  -  sin  £2  sin  w  cos  / 
rt ,  -  -  cos  £2  sin  t)  -  sii:  £2  cos  w  cos  / 
r„  -  sin  fit  cos  o 

r:i « sin  n  cos  e>  +  cos  £2  sin  <y  cos  i 
r. :  “  —  sin  Cl  sin  «  +  cos  £2  cos  c o  cos  t 
r:i «  -  cos  ft  sin  / 
r,, « sin  w  sin  / 

/  j; «  cos  w  sin  / 
r.v  *■  cos  / 

/  *  f||/|  *r  /',;(}  T  TijH 

J  m  r:ip  -u  ri:0  -t  r.j/T" 

A' 

2.  IJK  to  PQW  (inverse  of=l) 

P  -  r, ,7  +  r„J  +  riXK 
(J  *  /’|,/  "f  +  r,2A 

II  ■  /|j/  +  /jjA  +  /"jjA 

3.  IJK  to  RSW 

r„  *  cos  £2  cos  Ug  -  sin  £2  sin  u,  cos  i 
rx j « sin  £2  cos  Ug  +  sin  Ug  cos  £2  cos  / 
r,} «  sin  /  sin  «/„ 

«  -  cos  £2  sin  Ug  -  sin  £2  cos  «0  cos  / 

*  -  sin  £2  sin  i/0  +  cos  £2  cos  Ug  cos  / 
r:i «  cos  i/j  sin  i 
r„  =*  sin  £2  sin  / 
r};  =  -  cos  £2  sin  i 

l-ys  =  COj>  i 

t  t  /'.jA 

6  ~  T  /*../  T  NjA 

II  =  #*3j/  +  r3J  •}*  /  jjA 


SU 


4.  RSW  to  1JK  (inverse  of  =.') 

/  *  r„R  -r  r;t5  -r  r,.M‘ 

./ *r, ./{•»•  r..N  I  r 

A  *  t  rnff 

5.  !>QW  to  RSW 

r, |  ■  cos  v8 
r,j « sin  v0 
rn-U.U 
>•;.  »  -  sin  vs 

/':j  «  COS  V, 
r.j  ■  0.0 
/•„  -  0.0 
fj)  ■  0.0 
/„■  1.0 

lUrJ  +  ruQ  +  rn  I? 

r;i/^-r  r::Q  +  r:llT_ 

I*  “  >*||P  +  /}j(2  *r  *' ij lT* 

6.  RSW  to  PQW  (inverse  or  #5) 

/”  *■  ruR  *r  r:,S  + 

Q  -  >•,./(  -r  c;i  +  r,.ir 
ir-fuK  +  Hal  +  rair 


APPENDIX  C.  ORBITAL  ELEMENTS 


The  user  is  assumed  lo  be  studying  orbital  mechanics  and  should  understand  the 
orbital  elements  and  how  to  calculate  them.  A  brief  description  of  the  elements  and  the 
equations  used  to  calculate  the  elements  follow.  Por  a  detailed  explanation  of  the  clc- 
ments  and  the  equations  to  calculate  them  refer  to  Chapters  1  and  2  of  reference  1. 
Figure  6  on  page  S3  shows  the  orbital  elements  in  the  Geocentric-Equatorial  and 
perifocal  coordinate  system. 

1.  Angular  Momentum  ih): 

The  specific  ar.gular  momentum  is  a  constant  of  the  motion  of  the  satellite, 
defined  as  h  -  Fat. 

h  *  Fat  «  hjf  +  hjJ  +  It^K 

/'< "  >jvk  ~  Wj 
hj-W-Wk 
Ih  »  r.v,  -  r,v; 

h  °  \  hi  +  hj  +  K 

2.  Node  Vector  (n): 

The  node  vector  is  a  vector  pointing  along  the  line  of  nodes  in  the  direction 
of  the  ascending  node. 

n  —  Kxh  *  —hjl  -f  lui 

^  hj  -f  lif 

3.  Scmi-latus  rectum  (p): 

The  scmi-latus  rectum  is  a  geometric  constant  of  the  conic  section. 


A.  Eccentricity  ten 

The  eccentricity  is  a  constant  defining  the  shape  of  the  conic  orbit. 

F  =  tt  [O'2  -  t*  jF  -  (FF)F] 

e  = !  F  j 

5.  Semi-major  axis  (a): 


S2 


Figure  6.  Orbital  elements 


The  semi-major  axis  is  a  constant  defining  the  size  of  the  orbit. 
(l-t’:) 


6.  Inclination  (i): 

The  inclination  is  the  angle  between  the  'K'  unit  vector  in  the  'UK'  system  and 
the  angular  momentum  vector.  'h\ 


S3 


r  It-, 

i  ■  cos  H  — j— )  *  cos  1  — ) 

.  Longitude  of  ascending  node  (fir. 

The  longitude  ol‘  the  ascending  node  is  the  angle  in  the  fundamental  plane  , 
between  the  'I'  unit  vector  and  the  point  where  the  satellite  crosses  through  the 
fundamental  plane  in  a  northerly  direction  (ascending  node)  measured  counter¬ 
clockwise  when  viewed  from  the  north  side  of  the  fundamental  plane. 

n-cos  *(- f) 


S.  Argument  of  perigee  («): 

The  argument  of  perigee  is  the  angle  in  the  plane  of  the  satellite’s  orbit,  be¬ 
tween  the  ascending  node  and  the  perigee  point,  measured  in  the  direction  of  the 
satellite’s  motion. 


t?> 


w-'iitr) 


cos 


-l 


(lift  +  «;<?/) 

lie 


9.  True  anomaly  at  epoch  (v„): 

The  true  anomaly  at  epoch  is  the  angle  in  the  phne  of  the  satellite’s  orbit,  be¬ 
tween  perigee  and  the  position  of  the  satellite  at  a  particular  time.  i{.  called  the 
“epoch". 


Vo-cos"l(-^-) 

I  (J.  Argument  of  latitude  («,.): 

The  argument  of  latitude  is  the  angle  in  the  plane  of  the  orbit,  between  the 
ascending  node  and  the  radius  vector  to  the  satellite  at  time  /0. 


1 1.  Longitude  of  perigee  (11): 

The  longitude  of  perigee  is  the  angle  from  ’1'  to  perigee  measured  eastward  to 
the  ascending  node  and  then  in  the  orbital  plane  to  perigee. 

H  -  £1  ~  to 

12.  True  longitude  at  epoch  {/„): 

The  true  longitude  at  epoch  is  the  angle  between  T  and  r0  (the  radius  vector 
to  the  satellite  at  /„  measured  eastward  to  the  ascending  node  and  then  it.  the  orbital 
plane  to  r0. 


13.  Period  (per): 

The  period  is  the  time  the  for  the  satellite  to  complete  one  orbit.' 


Per  = 


•> 


P 


1-1.  Eccentric  anomaly  (EAl: 
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The  cceciuric  anomaly  is  the  ancle  between  the  perigee  and  a  position  on  an 
auxiliary  circle  circumscribed  about  the  ellipse  where  a  perpendicular  line  to  the 
major  axis  has  been  extended  from  the  epoch  location  of  the  satellite  to  the  auxil¬ 
iary  circle. 


EA  -  cos-1 


£  *r  cos(v) 

1  -r  e  COSlv) 


15.  Mean  motion  (n’l: 

The  mean  motion  is  defined  below: 


1(5.  Mean  anomaly  (MA): 

The  mean  anomaly  is  defined  below: 

MA  ■  //'(/  -  7)  *■  EA  —  c  sinfE.ll 


17.  Time  of  flight  (TF): 

The  time  of  flight  is  the  elapsed  time  from  when  the  satellite  was  at  perigee  to 
the  current  epoch. 


APPENDIX  D.  SAMPLE  ORBITS 


To  demonstrate  the  capabilities  of  the  program,  a  variety  of  orbital  plots  will  follow: 
L  Low  earth  orbit  (LEO). _ 

perifocal  coofttxxATz  system 


l#M4  J 

i «  4a.ooo  Dec.  a  -  mt  kmc-  0.100  pen  -  1.70  hours 

AVERAGE  KATZ  Of  CHAJ4CI  OP  ELEMENTS  PCX  SECOND 

DtArr  -  0.00  da/dt  -  0.00  dvdt-  0.00 

DMM/DT-  0.00  DMA/DT  ■  0.00  DLAN'DT-  0.00 
DH/DT-  0.00  DAJytrr-  0.00 

AVERAGE  MAGNITUDE  OP  PONCE*  PEX  UNIT  MASS  OCM/T**) 
CANTU  -  0.0  MOON  -  0 JO  SUN  «  0.0  DKAG  -  0.0 

PENIOEE  AfWfH 

RADIUS  (KM)  XP-  0900.0  KA  -  7044.4 

VELOCITY  OCM/SEC)  VP  -  SCI  VA  -  8.72 


Figure  7.  Unperturbed  Low  Earth  Orbit  (LEO) 


Figure  7  shows  the  perifocal  plot  of  a  satellite  in  an  unperturbed  low  earth 
orbit  (LEO).  The  initial  parameters  of  the  orbit  were  entered  as  follows: 
radius  of  perigee  (RP)  =  6500  km 
eccentricity  (c)  =  0.1 
inclination  (i)  =  45  degrees. 
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PERlfOCAL  COORDINATE  IWtM 


I  -  44.906  DEC.  A  -  7203.3  KM  E  -  0.066  PER  -  1.66  HOURS 

AVERAGE  RATE  Of  CHARGE  Of  ELEMENTS  PER  SECOND 
DI/DT-  4JK>*10'*  DX/DT-  6.36*1  O'*  DC/DT  ■  1.00*10"* 

DMM/DT  -  1.60*10'*  DMA/DY  •  SRIMO**  DLAH^DT "  6.46*10'f 
DH/DT-  0.63*10'*  DART7T  -  0.21*10'* 

AVERAGE  MAGNITUDE  Of  PORCES  fER  UNIT  MASS  (KM/TT) 
EARTH  -  0.6*10'*  MOON  -  0A*IO"*  *UK  -  4J*10*“  DRAG  -  1.4*10"* 

ft  t;  1C  EE  A po*** 

RADIUS  OCM)  Rf  -  6460.6  RA  -  7607 JO 

VELOCITY  (KM/5EO  Vf  -  6.20  VA  -  6.74 


Figure  8.  Perturbed  Low  Earth  Orbit  (LEO) 


With  perturbing  forces  applied  to  the  previous  LEO,  the  drag  force  will  be  the 
dominate  perturbing  force.  The  drag  will  act  as  a  negative-  velocity  change  applied 
in  the  area  of  perigee,  with  the  result  of  decreasing  the  semi-major  axis  length,  this 
in  effect  will  decrease  the  eccentricity  of  the  orbit,  as  can  be  seen  by  comparing  the 
orbital  data  of  the  unperturbed  LEO  in  Figure  7  on  page  S6  with  the  orbital  data 
of  the  perturbed  LEO  in  Figure  S. 
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2.  Circular  orbit. 


i 


1  -  30.007  DEC.  A  -  0000.1  KMC-  0.000  PER  -  1.01  HOURS 

AVERAGE  RATE  OF  CHANCE  OT  ELEMENTS  PER  SECOND 
Dl/DT-  4.02‘»0*’  DA/DT  »  0.74*10"  DE/DT-  3.03*10" 

DMM/DT  2-28*10"  DMA/DT  -  1-32*10"  DLAN/DT  *  7.32*10'’ 

DH/DT  -  5.00*10"  DAfyUT  ■  1J2*10"  « 

AVERACE  MAGNITUDE  OP  FORCES  PER  UNIT  MASS  (KM/***2) 

EARTH  "  1^*10"  MOON  «  0.0*1O'"  SUN  -  4.3*10*'*  DRAG  -  3.0*10"* 

PERIGEE  ApogM 

RADIUS  (KM)  RP  -  0004.4  RA  -  0007.0 

VELOCITY  (KM/SEC)  VP  -  7 JO  VA  -  7.53 

Figure  9.  Circular  Orbit 


An  example  of  the  plot  of  the  ground  track  of  a  sequence  of  three  60  degree 
inclined  perturbed  circular  orbits  with  a  radius  of  7000  km  is  shown  in  Figure  9. 
The  sequence  of  orbits  displays  the  precession  of  the  orbit  around  the  earth" 


3.  Transfer  orbit. 


ptRirocAL  coordinate  irrroi 


ttMIJ- 


n 


xw 


■M1J- 

1  »  C.000  OCC.  A  «  14310.9  KM  E  -  0.319  PER  -  4.93  HOURS 


AVER  ACE  RATE  OF  CHANCE  OF  ELEMENTS  PER  SECOND 
Di/mr-  o.oo  da/dt*  o.oo  wywr-  o.oo 
DMM/DT-  0.00  DMA/rrr  -  0.00  DLAHTrr  -  0.00 
DH/irr  -  o.oo  DAiytJT  -  o.oo 

AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KMA"2> 
EARTH  -  0.0  MOON  -  0.0  SUN  »  0.0  DRAG  -  0.0 

PERIGEE  Apof** 

RADIUS  OCM)  RP  -  7000.0  RA  -  22t»lS 

VELOCITY  (KM/SEC)  VP  -  0.30  VA  -  2: >3 


Figure  10.  Transfer  Orbit 


The  transfer  orbit  between  a  circular,  equatorial  LEO  and  a  molniya  orbit 
(high  eccentric  orbit)  is  shown  in  Figure  10.  A  velocity  increase  of  1.75  krn,'s  was 
applied  at  the  perigee  to  simulate  a  perigee  kick  to  boost  the  satellite  into  the 
molniya  orbit.  A  similar  velocity  change  could  then  be  applied  at  apogee  to  create 
a  high  altitude  circular  orbit,  or  a  n.gative  velocity  change  applied  at  perigee  could 
be  used  to  bring  the  satellite  back  to  a  LEO. 
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A.  Geosynchronous  orbit 


GROUND  TRACK 


1  *  #0.000  PEG.  A  *  4??M.6  KMC-  0.000  PER  -  23.00  HOURS 

AVERAGE  RATS' OT  CHANCE  Of  ELEMENTS  PER  SECOND 
DI/DT  *  I.W’ID''  DA/DT-  1.2«*10“*  DC/DT-  I.OO’IO'* 

DMM/DT-  3.26*10*"  DMA/DT  -  3.02*10’*  DLAN/DT  -  1.00*10’* 
DH/DT-  3.30*10“'  DAfVDT «  3.02*10'* 

AVERAGE  MAGNITUDE  OT  FORCES  PER  UNIT  MASS  (KM/S**2) 
EARTH  «  0.3*10'*  MOON  -  4-9*10’*  SUN  *  24J*I0’*  DRAG  -  S-riO’" 

PERIGEE  Apof«« 

RADIUS  (KM)  RP  -  422339  RA  -  42233.3 

VELOCITY  (k'KAEC)  VP  -  3.07  VA  -  3.07 


Figure  11.  Geosynchronous  Orbit 


The  ground  track  of  a  perturbed  geosynchronous  orbit  inclined  60  decrees  is 
shown  in  Figure  11.  The  orbit  displays  the  figure  eight  typical  with  inclined 
geosynchronous  orbits. 
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